diff --git a/.typos.toml b/.typos.toml index ad3ced22..55205726 100644 --- a/.typos.toml +++ b/.typos.toml @@ -9,4 +9,6 @@ ists = "ists" ispcs = "ispcs" eqs = "eqs" rhs = "rhs" -MTK = "MTK" \ No newline at end of file +MTK = "MTK" +# `Ot` is the transpose of the operator matrix `O` in operator inference. +Ot = "Ot" diff --git a/docs/src/functions.md b/docs/src/functions.md index 1fec4c49..89b24204 100644 --- a/docs/src/functions.md +++ b/docs/src/functions.md @@ -9,4 +9,10 @@ TSVD RSVD deim pod +baltrunc +BalancedTruncation +opinf +OperatorInferenceModel +reduced_dynamics +quadratic_monomials ``` diff --git a/docs/src/interfaces.md b/docs/src/interfaces.md index 7a1cf1d6..0a103cf3 100644 --- a/docs/src/interfaces.md +++ b/docs/src/interfaces.md @@ -2,7 +2,9 @@ The supported user workflow is deliberately small: construct a [`POD`](@ref) from state snapshots, then call [`reduce!`](@ref) with one of the built-in backends, -[`SVD`](@ref), [`TSVD`](@ref), or [`RSVD`](@ref). The public API does not require users +[`SVD`](@ref), [`TSVD`](@ref), or [`RSVD`](@ref). Continuous-time LTI systems can be +reduced with [`baltrunc`](@ref), and polynomial reduced models can be learned from +snapshot/derivative data with [`opinf`](@ref). The public API does not require users to interact with the abstract types described below. ## POD reduction contract diff --git a/src/ModelOrderReduction.jl b/src/ModelOrderReduction.jl index da03ed34..12e82c8b 100644 --- a/src/ModelOrderReduction.jl +++ b/src/ModelOrderReduction.jl @@ -3,11 +3,13 @@ module ModelOrderReduction using DocStringExtensions: DocStringExtensions, FUNCTIONNAME, TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES -using ModelingToolkit: ModelingToolkit, @parameters, @variables, Differential, Equation, Num, +using ModelingToolkit: ModelingToolkit, @parameters, @variables, Differential, Equation, + Num, System, SymbolicUtils, Symbolics, complete, mtkcompile, substitute using SciMLBase: SciMLBase using SymbolicIndexingInterface: SymbolicIndexingInterface -using LinearAlgebra: LinearAlgebra, /, \, mul!, qr, svd +using LinearAlgebra: LinearAlgebra, /, \, mul!, qr, svd, lyap, Symmetric, Diagonal, + eigen, eigvals, I, factorize, norm, ColumnNorm using SparseArrays: SparseArrays, SparseMatrixCSC, findnz, sparse include("Types.jl") @@ -24,6 +26,12 @@ export deim include("pod.jl") export pod +include("balanced_truncation.jl") +export BalancedTruncation, baltrunc + +include("operator_inference.jl") +export OperatorInferenceModel, opinf, reduced_dynamics, quadratic_monomials + include("precompile.jl") end diff --git a/src/balanced_truncation.jl b/src/balanced_truncation.jl new file mode 100644 index 00000000..7299057a --- /dev/null +++ b/src/balanced_truncation.jl @@ -0,0 +1,241 @@ +""" + $(TYPEDEF) + +Reduced continuous-time LTI model produced by [`baltrunc`](@ref). + +The reduced dynamics are +```math +\\dot x_r = A_r x_r + B_r u,\\qquad y = C_r x_r + D_r u, +``` +with reconstruction ``x \\approx T_r x_r`` and reduction ``x_r = S x`` for plain +truncation. When residualization is used, ``T_r`` is still the balanced truncation +map for the retained states; the static contribution of the discarded states is folded +into ``(A_r, B_r, C_r, D_r)`` rather than into ``T_r``. + +# Fields +$(TYPEDFIELDS) +""" +struct BalancedTruncation{T <: AbstractFloat} + "reduced state matrix ``A_r``" + A::Matrix{T} + "reduced input matrix ``B_r``" + B::Matrix{T} + "reduced output matrix ``C_r``" + C::Matrix{T} + "feedthrough matrix ``D_r``" + D::Matrix{T} + "retained Hankel singular values" + hsv::Vector{T} + "right projector ``T_r`` with ``x \\approx T_r x_r`` under plain truncation" + Tr::Matrix{T} + "left projector ``S`` with ``x_r = S x``" + S::Matrix{T} +end + +function Base.show(io::IO, bt::BalancedTruncation) + r = size(bt.A, 1) + nu = size(bt.B, 2) + ny = size(bt.C, 1) + return print(io, "BalancedTruncation(order = $r, inputs = $nu, outputs = $ny)") +end + +""" + $(TYPEDSIGNATURES) + +Return a square-root factor ``L`` of the positive-semidefinite matrix ``W`` such that +``W \\approx L L'``. Tiny negative eigenvalues from roundoff are clipped to zero. +""" +function _psd_sqrt(W::AbstractMatrix{T}; atol::Real = 0) where {T <: AbstractFloat} + F = eigen(Symmetric(Matrix{T}(W))) + λ = F.values + λmax = maximum(abs, λ) + thresh = max(T(atol), eps(T) * λmax) + kept = findall(≥(thresh), λ) + isempty(kept) && return zeros(T, size(W, 1), 0) + return F.vectors[:, kept] * Diagonal(sqrt.(λ[kept])) +end + +""" + $(TYPEDSIGNATURES) + +Compute continuous-time controllability and observability Gramians of ``(A,B,C)``. + +Solves +```math +A P_c + P_c A' + B B' = 0,\\qquad A' P_o + P_o A + C' C = 0 +``` +using `LinearAlgebra.lyap`. The open-loop matrix `A` must be Hurwitz. +""" +function _gramians(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) + Pc = lyap(A, B * B') + Po = lyap(A', C' * C) + return Symmetric(Pc), Symmetric(Po) +end + +function _select_order(σ::AbstractVector{T}, n, atol::Real, rtol::Real) where {T} + isempty(σ) && return 0 + if n !== nothing + n isa Integer || throw(ArgumentError("n must be an Integer or `nothing`")) + n < 0 && throw(ArgumentError("n must be nonnegative")) + return min(Int(n), length(σ)) + end + σmax = σ[1] + return count(s -> s >= T(atol) && s >= T(rtol) * σmax, σ) +end + +function _assert_hurwitz(A::AbstractMatrix{T}) where {T <: AbstractFloat} + if maximum(real, eigvals(A)) >= zero(T) + throw(ArgumentError("A must be Hurwitz (all eigenvalues with negative real part) for continuous-time balanced truncation")) + end + return nothing +end + +""" + baltrunc(A, B, C, D = nothing; n = nothing, atol = 0, rtol = 1e-3, + residual = false) -> BalancedTruncation + +Balanced truncation of a continuous-time linear time-invariant system +```math +\\dot x = A x + B u,\\qquad y = C x + D u. +``` + +The square-root method balances the controllability and observability Gramians and +truncates states with small Hankel singular values. If `n` is omitted, the reduced +order keeps singular values at least `atol` and at least `rtol` times the largest +singular value. + +When `residual = true`, truncated states are eliminated by residualization +(singular perturbation) so that the static gain is matched. + +Only real floating-point systems are supported. `A` must be Hurwitz. + +# Arguments +- `A::AbstractMatrix`: state matrix. +- `B::AbstractVecOrMat`: input matrix (vectors are treated as single-input). +- `C::AbstractMatrix`: output matrix. +- `D`: feedthrough matrix. Defaults to a zero matrix conforming to `(C, B)`. + +# Keywords +- `n = nothing`: reduced order. When `nothing`, choose the order from `atol` / `rtol`. +- `atol = 0`: absolute Hankel-singular-value cutoff. +- `rtol = 1e-3`: relative Hankel-singular-value cutoff. +- `residual = false`: use residualization instead of plain truncation. + +# Returns +- [`BalancedTruncation`](@ref): reduced matrices, retained Hankel singular values, and + projectors. + +# Throws +- `ArgumentError`: if the system dimensions are inconsistent, `A` is not square, `A` is + not Hurwitz, no positive Hankel singular values remain, or `n` is invalid. + +# Examples +```jldoctest +julia> using ModelOrderReduction + +julia> A = [-1.0 0.0; 0.0 -2.0]; B = reshape([1.0, 1.0], 2, 1); C = [1.0 1.0]; + +julia> bt = baltrunc(A, B, C; n = 1); + +julia> size(bt.A) +(1, 1) +``` + +# References +- Moore, B. C. (1981). Principal component analysis in linear systems. *IEEE + Transactions on Automatic Control*. +- Antoulas, A. C. (2005). *Approximation of Large-Scale Dynamical Systems*. SIAM. +""" +function baltrunc( + A::AbstractMatrix, + B::AbstractVecOrMat, + C::AbstractMatrix, + D::Union{AbstractMatrix, Nothing} = nothing; + n = nothing, + atol::Real = 0, + rtol::Real = 1.0e-3, + residual::Bool = false + ) + TA = float(promote_type(eltype(A), eltype(B), eltype(C))) + TA <: AbstractFloat || + throw(ArgumentError("baltrunc currently supports real floating-point systems only")) + A = Matrix{TA}(A) + B = Matrix{TA}(B isa AbstractVector ? reshape(B, length(B), 1) : B) + C = Matrix{TA}(C) + nx = size(A, 1) + size(A, 2) == nx || throw(ArgumentError("A must be square")) + size(B, 1) == nx || throw(ArgumentError("B must have the same number of rows as A")) + size(C, 2) == nx || throw(ArgumentError("C must have the same number of columns as A")) + if D === nothing + D = zeros(TA, size(C, 1), size(B, 2)) + else + D = Matrix{TA}(D) + size(D) == (size(C, 1), size(B, 2)) || + throw(ArgumentError("D must have size (size(C, 1), size(B, 2))")) + end + if n !== nothing + n isa Integer || throw(ArgumentError("n must be an Integer or `nothing`")) + n < 1 && throw(ArgumentError("n must be a positive Integer")) + end + + _assert_hurwitz(A) + Pc, Po = _gramians(A, B, C) + Lc = _psd_sqrt(Pc) + Lo = _psd_sqrt(Po) + (size(Lc, 2) == 0 || size(Lo, 2) == 0) && + throw(ArgumentError("controllability or observability Gramian has no positive eigenvalues; check stability and (A,B,C)")) + + U, σ, V = svd(Lo' * Lc; full = false) + # Drop numerically zero HSVs before forming balancing transforms. + σtol = max(eps(TA)^(2 / 3) * (isempty(σ) ? zero(TA) : σ[1]), TA(atol)) + keep = findall(>(σtol), σ) + isempty(keep) && + throw(ArgumentError("no Hankel singular values satisfy the truncation tolerances")) + U = U[:, keep] + σ = σ[keep] + V = V[:, keep] + + r = _select_order(σ, n, atol, rtol) + r == 0 && + throw(ArgumentError("no Hankel singular values satisfy the truncation tolerances")) + if n !== nothing && Int(n) > length(σ) + throw(ArgumentError("requested order n = $n exceeds the $(length(σ)) positive Hankel singular values")) + end + + i1 = 1:r + σr = σ[i1] + Σinvsqrt = Diagonal(inv.(sqrt.(σr))) + Tr = Lc * V[:, i1] * Σinvsqrt # x ≈ Tr xr + Sl = Σinvsqrt * U[:, i1]' * Lo' # xr = Sl x + + if residual && r < length(σ) + # Balanced realization on the retained positive-HSV subspace, then residualize. + Σinvsqrt_full = Diagonal(inv.(sqrt.(σ))) + Tfull = Lc * V * Σinvsqrt_full + Sfull = Σinvsqrt_full * U' * Lo' + Abal = Sfull * A * Tfull + Bbal = Sfull * B + Cbal = C * Tfull + i2 = (r + 1):length(σ) + A11 = Abal[i1, i1] + A12 = Abal[i1, i2] + A21 = Abal[i2, i1] + A22 = Abal[i2, i2] + B1 = Bbal[i1, :] + B2 = Bbal[i2, :] + C1 = Cbal[:, i1] + C2 = Cbal[:, i2] + A22fac = factorize(-A22) + A2221 = A22fac \ A21 + Ar = A11 + A12 * A2221 + Br = B1 + A12 * (A22fac \ B2) + Cr = C1 + C2 * A2221 + Dr = D + C2 * (A22fac \ B2) + return BalancedTruncation{TA}(Ar, Br, Cr, Dr, Vector{TA}(σr), Tr, Sl) + else + Ar = Sl * A * Tr + Br = Sl * B + Cr = C * Tr + return BalancedTruncation{TA}(Ar, Br, Cr, D, Vector{TA}(σr), Tr, Sl) + end +end diff --git a/src/operator_inference.jl b/src/operator_inference.jl new file mode 100644 index 00000000..c4e26f47 --- /dev/null +++ b/src/operator_inference.jl @@ -0,0 +1,390 @@ +""" + $(TYPEDEF) + +Reduced model learned by [`opinf`](@ref). + +With reduced coordinates ``x_r = V^\\top x``, the inferred dynamics are +```math +\\dot x_r = c + A x_r + B u + H\\,\\mathrm{quad}(x_r), +``` +where ``\\mathrm{quad}`` stacks the unique quadratic monomials +``x_i x_j`` for ``1 \\le i \\le j \\le r`` in row-major lower-triangular order +``(x_1^2, x_1 x_2, \\ldots, x_1 x_r, x_2^2, \\ldots, x_r^2)``. + +# Fields +$(TYPEDFIELDS) +""" +struct OperatorInferenceModel{T <: AbstractFloat} + "POD / trial basis ``V`` with orthonormal columns" + basis::Matrix{T} + "linear operator ``A``, or `nothing` if not inferred" + A::Union{Nothing, Matrix{T}} + "quadratic operator ``H`` with `size(H, 2) == r(r + 1) ÷ 2`, or `nothing`" + H::Union{Nothing, Matrix{T}} + "input operator ``B``, or `nothing` if not inferred" + B::Union{Nothing, Matrix{T}} + "constant term ``c``, or `nothing` if not inferred" + c::Union{Nothing, Vector{T}} +end + +function Base.show(io::IO, m::OperatorInferenceModel) + r = size(m.basis, 2) + terms = String[] + m.c === nothing || push!(terms, "c") + m.A === nothing || push!(terms, "A") + m.B === nothing || push!(terms, "B") + m.H === nothing || push!(terms, "H") + return print(io, "OperatorInferenceModel(order = $r, terms = $(join(terms, ", ")))") +end + +""" + $(TYPEDSIGNATURES) + +Fill `out` with the unique quadratic monomials of `x` in row-major lower-triangular +order ``(x_1^2, x_1 x_2, \\ldots, x_1 x_r, x_2^2, \\ldots, x_r^2)``. +""" +function quadratic_monomials!(out::AbstractVector{To}, x::AbstractVector{Tx}) where {To, Tx} + r = length(x) + length(out) == (r * (r + 1)) ÷ 2 || + throw(DimensionMismatch("quadratic monomial buffer has wrong length")) + k = 1 + @inbounds for i in 1:r + xi = x[i] + for j in i:r + out[k] = xi * x[j] + k += 1 + end + end + return out +end + +""" + $(TYPEDSIGNATURES) + +Allocate and return the unique quadratic monomials of `x`. +""" +function quadratic_monomials(x::AbstractVector{T}) where {T} + r = length(x) + out = Vector{float(T)}(undef, (r * (r + 1)) ÷ 2) + return quadratic_monomials!(out, x) +end + +""" + $(TYPEDSIGNATURES) + +Evaluate the reduced right-hand side of `model` at reduced state `xr` with optional +input `u`. +""" +function reduced_dynamics( + model::OperatorInferenceModel{T}, + xr::AbstractVector, + u::Union{AbstractVector, Nothing} = nothing + ) where {T} + length(xr) == size(model.basis, 2) || + throw(DimensionMismatch("reduced state has length $(length(xr)), expected $(size(model.basis, 2))")) + Tv = promote_type(T, eltype(xr)) + if u !== nothing + Tv = promote_type(Tv, eltype(u)) + end + dx = zeros(Tv, length(xr)) + if model.c !== nothing + dx .+= model.c + end + if model.A !== nothing + mul!(dx, model.A, xr, true, true) + end + if model.B !== nothing + u === nothing && + throw(ArgumentError("model has an input operator B but no input u was provided")) + length(u) == size(model.B, 2) || + throw(DimensionMismatch("input has length $(length(u)), expected $(size(model.B, 2))")) + mul!(dx, model.B, u, true, true) + end + if model.H !== nothing + q = quadratic_monomials(xr) + mul!(dx, model.H, q, true, true) + end + return dx +end + +""" + $(TYPEDSIGNATURES) + +Dense POD basis for Operator Inference. Prefer dense SVD over truncated SVD so that +requesting the full numerical rank still yields orthonormal columns. +""" +function _opinf_pod_basis(snapshot::AbstractMatrix{T}, dim::Integer) where {T} + reducer = POD(snapshot, dim) + reduce!(reducer, SVD()) + V = reducer.rbasis + if norm(V' * V - I(size(V, 2))) > sqrt(eps(T)) * size(V, 2) + throw(ArgumentError("computed POD basis is not orthonormal; try a smaller nmodes")) + end + return V +end + +function _assert_orthonormal_basis(V::AbstractMatrix{T}; atol = sqrt(eps(T))) where {T} + r = size(V, 2) + if norm(V' * V - I(r)) > atol * max(r, 1) + throw(ArgumentError("basis columns must be orthonormal")) + end + return nothing +end + +""" + $(TYPEDSIGNATURES) + +Build the Operator Inference data matrix whose columns are the selected operator +features for each snapshot. Snapshot data are stored with time along dimension 2 +in `Xr` (`r × k`). +""" +function _opinf_data_matrix( + Xr::AbstractMatrix{T}; + linear::Bool, + quadratic::Bool, + inputs::Union{Nothing, AbstractMatrix{<:Number}}, + constant::Bool + ) where {T} + r, k = size(Xr) + blocks = Matrix{T}[] + widths = Int[] + if constant + push!(blocks, ones(T, k, 1)) + push!(widths, 1) + end + if linear + push!(blocks, Matrix{T}(Xr')) + push!(widths, r) + end + if inputs !== nothing + U = Matrix{T}(inputs) + size(U, 2) == k || throw(ArgumentError("inputs must have one column per snapshot")) + push!(blocks, Matrix{T}(U')) + push!(widths, size(U, 1)) + end + if quadratic + nquad = (r * (r + 1)) ÷ 2 + Q = Matrix{T}(undef, k, nquad) + x = Vector{T}(undef, r) + q = Vector{T}(undef, nquad) + for j in 1:k + @inbounds for i in 1:r + x[i] = Xr[i, j] + end + quadratic_monomials!(q, x) + @inbounds for i in 1:nquad + Q[j, i] = q[i] + end + end + push!(blocks, Q) + push!(widths, nquad) + end + isempty(blocks) && + throw(ArgumentError("at least one of linear, quadratic, inputs, or constant must be enabled")) + return hcat(blocks...), widths +end + +function _unpack_opinf_operators( + O::AbstractMatrix{T}, + widths::Vector{Int}; + linear::Bool, + quadratic::Bool, + has_inputs::Bool, + constant::Bool + ) where {T} + offset = 0 + idx = 1 + c = nothing + A = nothing + B = nothing + H = nothing + if constant + w = widths[idx] + idx += 1 + c = vec(O[:, (offset + 1):(offset + w)]) + offset += w + end + if linear + w = widths[idx] + idx += 1 + A = O[:, (offset + 1):(offset + w)] + offset += w + end + if has_inputs + w = widths[idx] + idx += 1 + B = O[:, (offset + 1):(offset + w)] + offset += w + end + if quadratic + w = widths[idx] + H = O[:, (offset + 1):(offset + w)] + offset += w + end + offset == size(O, 2) || throw(ArgumentError("internal operator packing mismatch")) + return A, H, B, c +end + +function _normalize_inputs(inputs, k::Int, ::Type{T}) where {T} + inputs === nothing && return nothing + if inputs isa AbstractVector + length(inputs) == k || + throw(ArgumentError("vector inputs must have one entry per snapshot")) + return reshape(Vector{T}(inputs), 1, k) + end + U = Matrix{T}(inputs) + size(U, 2) == k || throw(ArgumentError("inputs must have one column per snapshot")) + return U +end + +""" + opinf(X, Xdot; nmodes, basis = nothing, linear = true, quadratic = false, + inputs = nothing, constant = false, λ = 0) -> OperatorInferenceModel + +Learn a polynomial reduced model by Operator Inference. + +Snapshot matrix `X` and derivative matrix `Xdot` store one snapshot per column +(state dimension × number of snapshots). A POD basis with `nmodes` columns is +computed from `X` unless `basis` is provided. The projected least-squares problem +```math +\\min_O \\| D O^\\top - \\dot X_r^\\top \\|_F^2 + \\lambda \\|O\\|_F^2 +``` +is solved for the selected operator blocks. + +# Arguments +- `X::AbstractMatrix`: state snapshots as columns. +- `Xdot::AbstractMatrix`: time-derivative snapshots as columns, conforming to `X`. + +# Keywords +- `nmodes::Integer`: number of POD modes. Required when `basis` is omitted. +- `basis = nothing`: optional orthonormal basis ``V``. When provided, `nmodes` defaults + to `size(basis, 2)`. +- `linear = true`: infer a linear operator ``A``. +- `quadratic = false`: infer a quadratic operator ``H`` on unique monomials. +- `inputs = nothing`: optional input snapshots as an `(input dimension × snapshots)` + matrix, or a length-`snapshots` vector for a single input. +- `constant = false`: infer a constant forcing term ``c``. +- `λ = 0`: Tikhonov regularization weight applied isotropically to all operator + coefficients. When positive, the ridge problem is solved via an augmented QR + factorization. + +# Returns +- [`OperatorInferenceModel`](@ref): inferred basis and operator blocks. + +# Throws +- `ArgumentError` / `DimensionMismatch`: if sizes are inconsistent, the basis is not + orthonormal, or no operator block is requested. + +# Examples +```jldoctest +julia> using ModelOrderReduction, LinearAlgebra + +julia> Atrue = [-1.0 0.0; 0.0 -2.0]; + +julia> X = reduce(hcat, [[exp(-t), exp(-2t)] for t in 0:0.05:3]); + +julia> Xdot = Atrue * X; + +julia> model = opinf(X, Xdot; basis = Matrix{Float64}(I, 2, 2)); + +julia> model.A ≈ Atrue +true +``` + +# References +- Peherstorfer, B. & Willcox, K. (2016). Data-driven operator inference for + nonintrusive projection-based model reduction. *CMAME*. +- Qian, E., Kramer, B., Peherstorfer, B. & Willcox, K. (2020). Lift & Learn. +""" +function opinf( + X::AbstractMatrix, + Xdot::AbstractMatrix; + nmodes::Union{Integer, Nothing} = nothing, + basis::Union{AbstractMatrix, Nothing} = nothing, + linear::Bool = true, + quadratic::Bool = false, + inputs = nothing, + constant::Bool = false, + λ::Real = 0 + ) + size(X) == size(Xdot) || throw(DimensionMismatch("X and Xdot must have the same size")) + T = float(promote_type(eltype(X), eltype(Xdot))) + T <: AbstractFloat || + throw(ArgumentError("opinf currently supports real floating-point snapshots only")) + X = Matrix{T}(X) + Xdot = Matrix{T}(Xdot) + nstate, k = size(X) + nstate >= 1 || throw(ArgumentError("state dimension must be positive")) + k >= 1 || throw(ArgumentError("need at least one snapshot")) + λ < 0 && throw(ArgumentError("regularization weight λ must be nonnegative")) + + if basis === nothing + nmodes isa Integer || + throw(ArgumentError("nmodes is required when basis is not provided")) + nmodes = Int(nmodes) + (0 < nmodes <= min(nstate, k)) || + throw(ArgumentError("nmodes must be in 1:min(state dimension, snapshots)")) + V = _opinf_pod_basis(X, nmodes) + else + V = Matrix{T}(basis) + size(V, 1) == nstate || throw(DimensionMismatch("basis must have $nstate rows")) + nmodes = nmodes === nothing ? size(V, 2) : Int(nmodes) + size(V, 2) == nmodes || + throw(ArgumentError("nmodes must equal the number of basis columns")) + nmodes >= 1 || throw(ArgumentError("basis must have at least one column")) + _assert_orthonormal_basis(V) + end + + Umat = _normalize_inputs(inputs, k, T) + Xr = V' * X + Rr = V' * Xdot + D, widths = _opinf_data_matrix( + Xr; + linear = linear, + quadratic = quadratic, + inputs = Umat, + constant = constant + ) + # D is k × p, Rr' is k × r; solve D * O' ≈ Rr' in the least-squares sense. + Rt = Matrix{T}(Rr') + p = size(D, 2) + if size(D, 1) < p && λ == 0 + @warn "Operator Inference least-squares problem is underdetermined (snapshots < features); consider more data, fewer terms, or λ > 0" snapshots = size( + D, 1 + ) features = p + end + if λ == 0 + Ot = qr(D, ColumnNorm()) \ Rt + else + √λ = sqrt(T(λ)) + Dλ = vcat(D, √λ * I(p)) + Rtλ = vcat(Rt, zeros(T, p, size(Rt, 2))) + Ot = qr(Dλ, ColumnNorm()) \ Rtλ + end + O = Matrix{T}(Ot') + A, H, B, c = _unpack_opinf_operators( + O, + widths; + linear = linear, + quadratic = quadratic, + has_inputs = Umat !== nothing, + constant = constant + ) + return OperatorInferenceModel{T}(V, A, H, B, c) +end + +""" + $(TYPEDSIGNATURES) + +Operator Inference from a vector-of-snapshots representation, matching the +[`POD`](@ref) snapshot convention. +""" +function opinf( + X::AbstractVector{<:AbstractVector}, + Xdot::AbstractVector{<:AbstractVector}; + kwargs... + ) + Xm = reduce(hcat, X) + Xdm = reduce(hcat, Xdot) + return opinf(Xm, Xdm; kwargs...) +end diff --git a/test/balanced_truncation.jl b/test/balanced_truncation.jl new file mode 100644 index 00000000..7568c347 --- /dev/null +++ b/test/balanced_truncation.jl @@ -0,0 +1,51 @@ +using Test +using LinearAlgebra +using ModelOrderReduction + +@testset "baltrunc: diagonal decay recovers dominant mode" begin + A = Diagonal([-1.0, -10.0, -100.0]) + B = reshape([1.0, 0.1, 0.01], 3, 1) + C = [1.0 0.1 0.01] + bt = baltrunc(Matrix(A), B, C; n = 1) + @test size(bt.A) == (1, 1) + @test length(bt.hsv) == 1 + @test bt.hsv[1] > 0 + @test size(bt.Tr) == (3, 1) + @test size(bt.S) == (1, 3) + @test bt.S * bt.Tr ≈ I(1) atol = 1.0e-10 + # Dominant pole near -1 + @test only(eigvals(bt.A)) ≈ -1.0 atol = 5.0e-2 +end + +@testset "baltrunc: order selection by atol/rtol" begin + A = Diagonal([-1.0, -2.0, -50.0]) + B = Matrix{Float64}(I, 3, 3) + C = Matrix{Float64}(I, 3, 3) + bt = baltrunc(Matrix(A), B, C; atol = 1.0e-8, rtol = 1.0e-2) + @test 1 <= size(bt.A, 1) <= 3 + @test issorted(bt.hsv; rev = true) +end + +@testset "baltrunc: residualization matches DC gain better than truncation" begin + A = [-1.0 0.0; 0.0 -10.0] + B = reshape([1.0, 1.0], 2, 1) + C = [1.0 1.0] + D = zeros(1, 1) + full_dc = only(-C * (A \ B) + D) + trunc = baltrunc(A, B, C, D; n = 1, residual = false) + resid = baltrunc(A, B, C, D; n = 1, residual = true) + trunc_dc = only(-trunc.C * (trunc.A \ trunc.B) + trunc.D) + resid_dc = only(-resid.C * (resid.A \ resid.B) + resid.D) + @test abs(resid_dc - full_dc) <= abs(trunc_dc - full_dc) + 1.0e-12 + @test resid_dc ≈ full_dc atol = 1.0e-10 +end + +@testset "baltrunc: dimension checks" begin + A = [-1.0 0.0; 0.0 -2.0] + B = reshape([1.0, 1.0], 2, 1) + C = [1.0 1.0] + @test_throws ArgumentError baltrunc(A, B, C; n = 0) + @test_throws ArgumentError baltrunc([1.0 2.0], B, C) + @test_throws ArgumentError baltrunc(A, reshape([1.0], 1, 1), C) + @test_throws ArgumentError baltrunc(Diagonal([0.5, -2.0]), B, C; n = 1) +end diff --git a/test/operator_inference.jl b/test/operator_inference.jl new file mode 100644 index 00000000..8c2e95ef --- /dev/null +++ b/test/operator_inference.jl @@ -0,0 +1,136 @@ +using Test +using LinearAlgebra +using ModelOrderReduction + +@testset "opinf: recovers linear diagonal dynamics" begin + Atrue = Diagonal([-1.0, -2.0, -3.0]) + ts = 0:0.05:4 + X = reduce(hcat, [[exp(Atrue[i, i] * t) for i in 1:3] for t in ts]) + Xdot = Matrix(Atrue) * X + model = opinf(X, Xdot; basis = Matrix{Float64}(I, 3, 3)) + @test size(model.basis) == (3, 3) + @test model.A ≈ Matrix(Atrue) atol = 1.0e-8 + @test model.H === nothing + @test model.B === nothing + @test model.c === nothing + xr = X[:, 1] + @test reduced_dynamics(model, xr) ≈ model.A * xr + + pod_model = opinf(X, Xdot; nmodes = 2) + @test size(pod_model.basis) == (3, 2) + @test size(pod_model.A) == (2, 2) + Xr = pod_model.basis' * X + Rr = pod_model.basis' * Xdot + @test norm(Rr - pod_model.A * Xr) / norm(Rr) < 1.0e-2 +end + +@testset "opinf: recovers quadratic system with identity basis" begin + # ẋ₁ = -2 x₁ + x₁ x₂ + # ẋ₂ = -x₂ + x₂² + Atrue = [-2.0 0.0; 0.0 -1.0] + # unique monomials [x1^2, x1*x2, x2^2] + Htrue = [0.0 1.0 0.0; 0.0 0.0 1.0] + rng_states = [ + Float64[x1, x2] + for x1 in range(-0.5, 0.5; length = 9), + x2 in range(-0.5, 0.5; length = 9) + ] + X = reduce(hcat, vec(rng_states)) + Xdot = similar(X) + for j in axes(X, 2) + x = X[:, j] + Xdot[:, j] = Atrue * x + Htrue * quadratic_monomials(x) + end + model = opinf( + X, Xdot; basis = Matrix{Float64}(I, 2, 2), linear = true, quadratic = true + ) + @test model.A ≈ Atrue atol = 1.0e-10 + @test model.H ≈ Htrue atol = 1.0e-10 + x = [0.2, -0.3] + @test reduced_dynamics(model, x) ≈ Atrue * x + Htrue * quadratic_monomials(x) +end + +@testset "opinf: forced linear system with inputs" begin + Atrue = [-1.0 0.0; 0.0 -2.0] + Btrue = reshape([1.0, 0.5], 2, 1) + ts = 0:0.02:2 + U = reshape([sin(4t) for t in ts], 1, :) + # Integrate exactly for diagonal A with forcing: use discrete update from known ODE solution + # Generate consistent (X, Xdot, U) pairs on a grid of states instead. + xs = [ + Float64[x1, x2] + for x1 in range(-1, 1; length = 7), x2 in range(-1, 1; length = 7) + ] + us = range(-1, 1; length = 5) + cols_x = Vector{Float64}[] + cols_dx = Vector{Float64}[] + cols_u = Vector{Float64}[] + for x in vec(xs), u in us + + push!(cols_x, x) + push!(cols_dx, Atrue * x + Btrue * [u]) + push!(cols_u, [u]) + end + X = reduce(hcat, cols_x) + Xdot = reduce(hcat, cols_dx) + Umat = reduce(hcat, cols_u) + model = opinf(X, Xdot; basis = Matrix{Float64}(I, 2, 2), inputs = Umat) + @test model.A ≈ Atrue atol = 1.0e-10 + @test model.B ≈ Btrue atol = 1.0e-10 +end + +@testset "opinf: vector-of-snapshots API and regularization" begin + Xcols = [[exp(-0.5t), exp(-1.5t)] for t in 0:0.1:2] + Xdotcols = [[-0.5 * exp(-0.5t), -1.5 * exp(-1.5t)] for t in 0:0.1:2] + model = opinf(Xcols, Xdotcols; nmodes = 2, λ = 1.0e-12) + @test size(model.A) == (2, 2) + @test size(model.basis, 2) == 2 +end + +@testset "opinf: constant term and input evaluation" begin + Atrue = [-1.0 0.0; 0.0 -2.0] + ctrue = [0.25, -0.5] + Btrue = reshape([1.0, 0.5], 2, 1) + xs = [ + Float64[x1, x2] + for x1 in range(-1, 1; length = 6), x2 in range(-1, 1; length = 6) + ] + us = range(-1, 1; length = 5) + cols_x = Vector{Float64}[] + cols_dx = Vector{Float64}[] + cols_u = Float64[] + for x in vec(xs), u in us + + push!(cols_x, x) + push!(cols_dx, Atrue * x + Btrue * [u] + ctrue) + push!(cols_u, u) + end + X = reduce(hcat, cols_x) + Xdot = reduce(hcat, cols_dx) + model = opinf( + X, + Xdot; + basis = Matrix{Float64}(I, 2, 2), + inputs = cols_u, + constant = true + ) + @test model.A ≈ Atrue atol = 1.0e-10 + @test model.B ≈ Btrue atol = 1.0e-10 + @test model.c ≈ ctrue atol = 1.0e-10 + x = [0.2, -0.1] + u = [0.3] + @test reduced_dynamics(model, x, u) ≈ Atrue * x + Btrue * u + ctrue + @test_throws ArgumentError reduced_dynamics(model, x) +end + +@testset "opinf: quadratic monomials for r=3" begin + x = [2.0, 3.0, 5.0] + @test quadratic_monomials(x) ≈ [4.0, 6.0, 10.0, 9.0, 15.0, 25.0] +end + +@testset "opinf: argument checks" begin + X = rand(3, 5) + Xdot = rand(3, 4) + @test_throws DimensionMismatch opinf(X, Xdot; nmodes = 2) + @test_throws ArgumentError opinf(rand(3, 5), rand(3, 5); nmodes = 2, linear = false) +end