B/polynomialize quadratize reduce

.gitignore

@@ -23,4 +23,8 @@ docs/site/
 # environment.
 Manifest.toml
 
-.vscode
+.vscode
+
+# macOS metadata
+.DS_Store
+**/.DS_Store

Project.toml

@@ -7,6 +7,7 @@ version = "0.1.6"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 RandomizedLinAlg = "0448d7d9-159c-5637-8537-fd72090fea46"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
@@ -17,6 +18,7 @@ TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e"
 DocStringExtensions = "0.8, 0.9"
 LinearAlgebra = "1"
 ModelingToolkit = "11"
+OrdinaryDiffEq = "6, 7"
 PrecompileTools = "1"
 RandomizedLinAlg = "0.1"
 Setfield = "0.8, 1"

README.md

@@ -81,3 +81,46 @@ sol_deim_w = deim_sol[w(x, t)]
 The following figure shows the comparison of the solutions of the 32-dimension full-order model and the POD5-DEIM5 reduced-order model.
 
 ![comparison](https://user-images.githubusercontent.com/45696147/195765614-df9092a2-4fca-4602-bb15-81e65b2b572e.svg)
+
+#### Polynomialization, Quadratization, and Galerkin Reduction on an ODE system
+```julia
+using ModelOrderReduction
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+using OrdinaryDiffEq
+
+#Create the ModelingToolkit System
+@variables x(t) y(t)
+
+eqs = [
+    D(x) ~ -x + y + 0.1 * sqrt(x),
+    D(y) ~ -2.0 * y + 0.2 * x^2,
+]
+
+@mtkcompile sys = System(eqs, t)
+
+#Provide initial conditions and time span for ODE
+u0 = [1.0, 0.5]
+tspan = (0.0, 1.0)
+
+#Provide number of variables in reduced ODE
+nmodes = 1
+
+result = polynomialize_quadratize_reduce(
+    sys,
+    u0,
+    tspan,
+    nmodes;
+    saveat = 0.1,
+)
+
+#result has a number of components, in particular:
+#rom: a ModelingToolkit System for the reduced order model
+@show result.rom
+#a0: initial conditions for the reduced model
+@show result.a0
+#V: basis used for the affine Galerkin reduction
+@show size(result.V)
+#xbar: center used for the affine Galerkin reduction
+@show size(result.xbar)
+```

src/ModelOrderReduction.jl

@@ -2,9 +2,10 @@ module ModelOrderReduction
 
 using DocStringExtensions: DocStringExtensions, FUNCTIONNAME, SIGNATURES, TYPEDSIGNATURES
 
-using ModelingToolkit: ModelingToolkit, @variables, Differential, Equation, Num, ODESystem,
-    SymbolicUtils, Symbolics, arguments, build_function, complete,
-    expand, substitute, tearing_substitution
+using ModelingToolkit: ModelingToolkit, @variables, @named, Differential, Equation, Num,
+    ODESystem, System, SymbolicUtils, Symbolics, arguments, build_function, complete,
+    equations, expand, iscall, operation, substitute, tearing_substitution, unknowns
+using OrdinaryDiffEq: ODEProblem, Tsit5, solve
 using LinearAlgebra: LinearAlgebra, /, \, mul!, qr, svd
 
 using Setfield: Setfield, @set!
@@ -21,6 +22,18 @@ export POD, reduce!
 include("deim.jl")
 export deim
 
+include("PolynomializeQuadratizeReduce/PolynomializeQuadratizeReduceUtils.jl")
+include("PolynomializeQuadratizeReduce/Polynomialization.jl")
+include("PolynomializeQuadratizeReduce/Quadratization.jl")
+include("PolynomializeQuadratizeReduce/GalerkinReduction.jl")
+include("PolynomializeQuadratizeReduce/PolynomializeQuadratizeReduce.jl")
+
+export polynomialize
+export quadratize
+export galerkin_project_system
+export galerkin_project_system_affine
+export polynomialize_quadratize_reduce
+
 include("precompile.jl")
 
 end

src/PolynomializeQuadratizeReduce/GalerkinReduction.jl

@@ -0,0 +1,155 @@
+function differentiated_variable(eq)
+    lhs = unwrap(eq.lhs)
+    args = arguments(lhs)
+    return Num(args[1])
+end
+
+function ordered_rhs(sys)
+    eqs = equations(sys)
+    xs = unknowns(sys)
+
+    rhs_by_var = Dict{Any, Any}()
+
+    for eq in eqs
+        x = differentiated_variable(eq)
+        rhs_by_var[x] = eq.rhs
+        rhs_by_var[unwrap(x)] = eq.rhs
+    end
+
+    return [rhs_by_var[x] for x in xs]
+end
+
+"""
+    galerkin_project_system_affine(sys, V, xbar, a_vars; pmap=Dict(), name=:rom)
+
+Project an explicit ODE system onto the affine trial space
+
+    x(t) ≈ xbar + V*a(t).
+
+Here `x(t)` is the full state vector of `sys`, `xbar` is a fixed offset vector,
+`V` is an `n × r` basis matrix, and `a(t)` is the reduced state vector with
+entries `a_vars`.
+
+The reduced system is constructed by direct Galerkin projection:
+
+    a'(t) = V' * f(xbar + V*a(t)),
+
+where `f` is the right-hand side of the full system. The returned system also
+contains observed equations reconstructing the full state variables as
+
+    x_i(t) ~ xbar[i] + sum(V[i, α] * a_α(t) for α in 1:r).
+
+Arguments:
+- `sys`: ModelingToolkit ODE system.
+- `V`: projection basis of size `(n, r)`, where `n = length(unknowns(sys))`.
+- `xbar`: affine offset vector of length `n`.
+- `a_vars`: reduced state variables of length `r`.
+
+Keywords:
+- `pmap`: optional parameter substitutions applied before projection.
+- `name`: name of the returned reduced system.
+
+Returns:
+- A ModelingToolkit `System` for the affine Galerkin reduced-order model.
+"""
+function galerkin_project_system_affine(sys, V, xbar, a_vars; pmap = Dict(), name = :rom)
+    xs = unknowns(sys)
+
+    n = length(xs)
+    n == size(V, 1) || error("V must have size (n,r), where n = length(unknowns(sys))")
+    length(xbar) == n || error("xbar must have length n")
+    r = size(V, 2)
+    length(a_vars) == r || error("length(a_vars) must equal size(V,2)")
+
+    iv = ModelingToolkit.get_iv(sys)
+    Dred = Differential(iv)
+
+    rhs = unwrap.(ordered_rhs(sys))
+
+    if !isempty(pmap)
+        subdict = Dict(Symbolics.unwrap(k) => v for (k, v) in pmap)
+        rhs = Symbolics.substitute.(rhs, Ref(subdict))
+    end
+
+    x_subs = Dict{Any, Any}()
+
+    for i in 1:n
+        rec = xbar[i]
+
+        for α in 1:r
+            rec += V[i, α] * a_vars[α]
+        end
+
+        x_subs[xs[i]] = rec
+        x_subs[unwrap(xs[i])] = rec
+    end
+
+    f_affine = [
+        Symbolics.substitute(rhs[i], x_subs)
+            for i in 1:n
+    ]
+
+    rhs_red = Vector{Any}(undef, r)
+
+    for α in 1:r
+        expr = zero(Num)
+
+        for i in 1:n
+            expr += V[i, α] * f_affine[i]
+        end
+
+        rhs_red[α] = Symbolics.simplify(Symbolics.expand(expr))
+    end
+
+    eqs_red = [Dred(a_vars[α]) ~ rhs_red[α] for α in 1:r]
+
+    obs = Vector{Equation}(undef, n)
+
+    for i in 1:n
+        rec = xbar[i]
+
+        for α in 1:r
+            rec += V[i, α] * a_vars[α]
+        end
+
+        obs[i] = xs[i] ~ Symbolics.simplify(rec)
+    end
+
+    return System(eqs_red, iv; observed = obs, name = name)
+end
+
+"""
+    galerkin_project_system(sys, V, a_vars; pmap=Dict(), name=:rom)
+
+Project an explicit ODE system onto the linear trial space
+
+    x(t) ≈ V*a(t).
+
+This is the zero-offset special case of `galerkin_project_system_affine`, namely
+`xbar = zeros(n)`, where `n = length(unknowns(sys))`.
+
+Arguments:
+- `sys`: ModelingToolkit ODE system.
+- `V`: projection basis of size `(n, r)`.
+- `a_vars`: reduced state variables of length `r`.
+
+Keywords:
+- `pmap`: optional parameter substitutions applied before projection.
+- `name`: name of the returned reduced system.
+
+Returns:
+- A ModelingToolkit `System` for the linear Galerkin reduced-order model.
+"""
+function galerkin_project_system(sys, V, a_vars; pmap = Dict(), name = :rom)
+    n = length(unknowns(sys))
+    xbar = zeros(Float64, n)
+
+    return galerkin_project_system_affine(
+        sys,
+        V,
+        xbar,
+        a_vars;
+        pmap = pmap,
+        name = name,
+    )
+end