GitHub

ModelOrderReductionToolkit.jl Models and Reductors

Stationary Models

A model <: StationaryModel{NOUT} is a map from parameter space to a vector space. When calling model(p, i=1) for i=1,...,NOUT, must return a vector of length output_length(model) and of eltype output_type(model).

Linear Model

A model of the form

\[A(p) x(p) = b(p).\]

ModelOrderReductionToolkit.LinearModel — Type
model = LinearModel(Ap::APArray, bp::APArray) <: StationaryModel{1}
model = LinearModel(Ap::APArray, b::AbstractVector) <: StationaryModel{1}

Struct for containing a parameterized linear model A(p) x = b(p) or A(p) x = b with affine parameter dependence. Can form a solution for a new parameter value by calling it on a new parameter value x = model(p).

source

For an example, see PoissonModel().

Linear Matrix Model

A model of the form

\[A(p) X(p) = B(p).\]

ModelOrderReductionToolkit.LinearMatrixModel — Type

LinearMatrixModel{NOUT}

model = LinearMatrixModel(Ap::APArray, bps::AbstractVector{Union{APArray, <:AbstractVector}})
model = LinearMatrixModel(Ap::APArray, Bp::APArray)
model = LinearMatrixModel(Ap::APArray, B::AbstractMatrix)

Struct for containing a parameterized linear model A(p) X = B(p), with affine parameter dependence. Stored internally similarly to a LinearModel, so can solve for a single column of the solution by model(p, i) for i=1,...,NOUT, or can solve for the matrix solution by model(p, 0).

source

For an example, see to_frequency_domain(PenzlModel()).

Stationary Model Reductors

Reduced basis methods for stationary models typically require the models to have implemented galerkin_project(model, V[, W=V], r=-1) and galerkin_add!(rom, fom, v, Vold[, w=v, Wold=Vold]).

POD Reductor

Proper orthogonal decomposition which solves the full order model at each parameter value, computes an SVD of the snapshot matrix, and allows for projection onto the resulting left singular vectors.

ModelOrderReductionToolkit.PODReductor — Type

pod_reductor <: PODReductor

A struct for holding the parts of a POD reductor for a StationaryModel model. Can access the FOM through pod_reductor.model, the snapshot matrix through pod_reductor.snapshots, the singular values through pod_reductor.S, and the reduced basis (left singular vectors) through pod_reductor.V.

source
model = PoissonModel()
params = [[i,j,k] for i in range(0,1,5) for j in range(0,1,5) for k in range(0,1,5)]
reductor = PODReductor(model)
add_to_rb!(reductor, params) # or add_to_rb!(reductor, snapshots) if snapshot matrix already formed
rom = form_rom(reductor, 10)

SG Reductor

Strong greedy method which solves the full order model at each parameter value, computes a column-pivoted QR decomposition of the snapshot matrix, and allows for projection onto the resulting orthonormalized snapshots.

ModelOrderReductionToolkit.SGReductor — Type

sg_reductor <: SGReductor

A struct for holding the parts of an SG reductor for a StationaryModel model. Can access the FOM through sg_reductor.model, the snapshot matrix through sg_reductor.snapshots, the pivot order through sg_reductor.p, and the reduced basis through sg_reductor.V.

source
model = PoissonModel()
params = [[i,j,k] for i in range(0,1,5) for j in range(0,1,5) for k in range(0,1,5)]
reductor = SGReductor(model)
add_to_rb!(reductor, params) # or add_to_rb!(reductor, snapshots) if snapshot matrix already formed
rom = form_rom(reductor, 10)

WG Reductor

Weak greedy method which utilizes an ErrorEstimator to attempt to choose the parameter value (and index) for which the error is the worst for the given RB. Then, forms that full order solution and continues.

ModelOrderReductionToolkit.WGReductor — Type

wg_reductor <: WGReductor

Stores a FOM in model, an error estimator estimator, the greedily selected parameters params_greedy, the reduced basis in V, the ROM in rom, the approximate errors at each step in approx_errors, and the truth errors in each step at truth_errors.

source
model = PoissonModel()
params = [[i,j,k] for i in range(0,1,5) for j in range(0,1,5) for k in range(0,1,5)]
scm = SCM(model.Ap, params, coercive=true)
estimator = StabilityResidualErrorEstimator(model, scm)
reductor = WGReductor(model, estimator)
add_to_rb!(reductor, params, 10)
rom = form_rom(reductor, 10)

See sources 1-3 on more information on successive constraint methods, residual computations, and methodology on weak greedy reduced basis methods.

Nonstationary Models

A model <: NonstationaryModel is a map from parameter space to a set of ODEs. Calling model(p) initializes the model to the given parameter value. ODE outputs are vector of length output_length(model) and of eltype output_type(model). Can convert into an ODEProblem by to_ode_problem(model[, x0, tspan, p]).

LTI Model

A model of the form

\[E(p) x'(t;p) = A(p) x(t;p) + B(p) u(t)\\ y(t;p) = C(p) x(t;p) + D(p) u(t)\]

ModelOrderReductionToolkit.LTIModel — Type

model = LTIModel(A_in, B_in, C_in, D_in=0, E_in=I) <: NonstationaryModel model = LTIModel(lti<:AbstractStateSpace) <: NonstationaryModel model = LTIModel(lti<:AbstractDescriptorStateSpace) <: NonstationaryModel

Struct for containing a parameterized LTI model E(p) x'(t,p) = A(p) x(t,p) + B(p) u(t) y(t,p) = C(p) x(t,p) + D(p) u(t) with affine parameter dependence where u(t) is some input signal. Can initialize the system to a given parameter p by calling model(p).

source

Also implements to_ss(model[, p=nothing]) to convert a model to a ControlSystems.jl state space, and to_dss(model[, p=nothing]) to convert to a DescriptorSystems.jl. Also implements to_frequency_domain(model) which returns a LinearMatrixModel to solve for the state variable x in the frequency domain. Also has galerkin_project(model, V[, W=V; r=-1]) implemented for forming reduced order models. To cast to an ODE problem for a given input u(x), call to_ode_problem(model[, x0=0.0, tspan=(0,1), p=nothing, u]).

BT Reductor

The state of the art method for reducing a non-parameterized LTI problem is through balanced truncation. This can be performed with a BTReductor object.

ModelOrderReductionToolkit.BTReductor — Type

reductor = BTReductor(model::LTIModel[, p=nothing; noise=0, iterative=nothing, maxdim=-1, lradi_eps=1e-6, dense_row_tol=1e-8])

Balanced truncation reductor object for reducing an LTIModel. If parameter p passed in, model first initialized to parameter value. noise determines amount of printed output. If isnothing(iterative), checks the size and sparsity of the system whether or not to use an iterative method. If iterative==true, uses an iterative low-rank ADI method to solve for the Gramians (see Kurschner and Benner 2016 Dissertation Alg 4.3), otherwise, uses MatrixEquations.jl to solve them densely. maxdim determines the maximum rank for the iterative solver, and lradi_eps determines a tolerance for when the algorithm terminantes. If iterative==false, dense_row_tol is used to truncate the dense Gramians for quicker solving of the RB and HSVs.

The reachability Gramian can be computed by reductor.R * reductor.R', and the observability Gramian by reductor.L * reductor.L'.

HSVs stored in reductor.hs.

Petrov-Galerkin test and trial spaces stored in reductor.V and reductor.W respectively.

Initialized-to parameter value stored in reductor.p.

source

After forming a BTReductor object on an LTIModel, can obtain the system Gramians through reachability_gramian(reductor) and observability_gramian(reductor). As with other reductors, has form_rom and lift implemented.

model = PenzlModel()
reductor = BTReductor(model)
rom = form_rom(reductor, 20)

See source 4 for more information on the iterative method for solving Lyapunov equations used by BTReductor when iterative==true, and see MatrixEquations.jl for the non-sparse Lyapunov solver. See source 5 for the Penzl model example and for more information on truncation of LTI systems.

RB Reduction

For reducing a parameterized LTI problem in a reduced basis setting, one option is to create a (complex) basis in the frequency domain.

model = ParameterizedPenzlModel()
freq_model = to_frequency_domain(model)
reductor = PODReductor(freq_model)
params = [[ω,p1,p2,p3] for ω in range(-500,500,21) for p1 in range(-25,25,6) for p2 in range(-25,25,6) for p3 in range(-25,25,6)]
reductor = PODReductor(freq_model)
add_to_rb!(reductor, params)
rom = galerkin_project(model, Matrix(reductor.V[:,1:20])) # Faster when converting from VOV to Matrix

References:

  1. D.B.P. Huynh, G. Rozza, S. Sen, A.T. Patera. A successive constraint linear optimization method for lower bounds of parametric coercivity and inf–sup stability constants. Comptes Rendus Mathematique. Volume 345, Issue 8. 2007. Pages 473-478. https://doi.org/10.1016/j.crma.2007.09.019.
  2. Quarteroni, Alfio, Andrea Manzoni, and Federico Negri. Reduced Basis Methods for Partial Differential Equations. Vol. 92. UNITEXT. Cham: Springer International Publishing, 2016. http://link.springer.com/10.1007/978-3-319-15431-2.
  3. Yanlai Chen, Jiang Jiahua, and Akil Narayan. A robust error estimator and a residual-free error indicator for reduced basis methods. Computers & Mathematics with Applications. 2019. http://www.sciencedirect.com/science/article/pii/S0898122118306850
  4. Patrick Kürschner and Peter Benner. Efficient low-rank solutions of large-scale matrix equations. Forschungsberichte aus dem Max-Planck-Institut für Dynamik Komplexer Technischer Systeme. 2016. https://pure.mpg.de/rest/items/item22467967/component/file_2296741/content.
  5. Thilo Penzl. Algorithms for model reduction of large dynamical systems. Linear Algebra and its Applications. Volume 415, Issue 2, Pages 322-343. June 1, 2006. https://www.sciencedirect.com/science/article/pii/S0024379506000371.
  6. D.B.P. Huynh, D.J. Knezevic, Y. Chen, J.S. Hesthaven, A.T. Patera. A natural-norm Successive Constraint Method for inf-sup lower bounds. Computer Methods in Applied Mechanics and Engineering. Volume 199, Issue 29. June 1, 2010. https://www.sciencedirect.com/science/article/pii/S0045782510000691.