Additional ModelOrderReductionToolkit.jl Docstrings

rom = galerkin_project(model::LinearModel, V[, W=V; r=-1])

Perform (Petrov) Galerkin projection on a linear model where the trial space is the first r columns of V and the test space is the first r columns of W. If r=-1, uses all columns. Returns a new LinearModel.

W' * A(p) * V * x_r = W' * b(p),V * x_r ≈ x = A(p)^(-1) b(p)`

rom = galerkin_project(model::LinearMatrixModel, V[, W=V; r=-1])

Perform (Petrov) Galerkin projection on each linear model in model.models.

galerkin_project(model, V[, W=V; WTEVisI=false, r=-1])

Performs Galerkin projection on the model <: LTIModel and returns a new LTIModel. By default, assumes that WᵀV=I, if WTEVisI==true, then assumes that WᵀEV=I.

galerkin_add!(rom::LinearModel, fom::LinearModel, v, Vold[, w=v, Wold=Vold; r=-1])

Assuming that rom = galerkin_project(model, Vold, Wold), updates rom such that if V = [Vold v] and W = [Wold w], then rom = galerkin_project(model, V, W).

galerkin_add!(rom::LinearMatrixModel, fom::LinearMatrixModel, v, Vold[, w=v, Wold=Vold; r=-1])

Assuming that rom = galerkin_project(fom, Vold, Wold), updates rom such that if V = [Vold v] and W = [Wold w], then rom = galerkin_project(fom, V, W).

model = PoissonModel([; Nx=999, P=3])

Uses finite differences on the following PDE to generate a LinearModel with parameter dependence. Change discretization with Nx, and change number of parameters, and hence number of affine terms, with P.

- ∂_x (κ(x, p) ∂_x u(x, p)) = f(x)

u(0,p) = 0; u(1,p) = p[1]

κ(x,p) = (1.05 - (1/2)^(P-1)) + 0.5 p[2] sin(2πx) + 0.25 p[3] sin(4πx) + ... + (1/2)^(P-1) p[P] sin(2π * P * x)

f(x,p) = (0.25 .< x .< 0.75) .* 10.0

length(p) = P; 0 ≤ p[i] ≤ 1

model = PenzlModel([ns=1006])

Generates the standard Penzl LTIModel with one input, one output, and ns=1006 dimension state variable.

model = MISOPenzlModel([ns=1006])

Generates an LTIModel with three inputs, one output, and a state of dimension ns=1006. Same structure as the Penzl model except the B matrix is changed to

1006×3 Matrix{Float64}:
 10.0   0.0   0.0
 10.0   0.0   0.0
  0.0  10.0   0.0
  0.0  10.0   0.0
  0.0   0.0  10.0
  0.0   0.0  10.0
  1.0   1.0   1.0
  ⋮
  1.0   1.0   1.0
  1.0   1.0   1.0
  1.0   1.0   1.0
  1.0   1.0   1.0
  1.0   1.0   1.0
  1.0   1.0   1.0
  1.0   1.0   1.0

model = ParameterizedPenzlModel([ns=1006])

Generates an LTIModel with one input, one output, and a state of dimension ns=1006. Same structure as the Penzl model, expect depends on a parameter vector of length 3 which shift the poles along the complex axis. Instantiate to a parameter vector by calling model([p1,p2,p3]).

frequency_model = to_frequency_domain(model, logfreq=false)

Assuming null initial conditions, uses the Laplace transform to convert the LTIModel into a LinearMatrixModel for which the first element of the parameter vector is the imaginary part of the frequency variable. The model is of the form sE(p) X = A(p) X + B(p) where X is the Laplace variable for X and s = 0 + iω. If logfreq=true, then scales ω such that p[1] = log10(ω).

to_ss(model[, p=nothing])

Must import ControlSystems.jl for this functionality as its definition lives in an extension.

Initializes the model to the parameter p if passed in, then returns a ControlSystems.jl StateSpace object.

to_dss(model[, p=nothing])

Must import DescriptorSystems.jl for this functionality as its definition lives in an extension.

Initializes the model to the parameter p if passed in, then returns a DescriptorSystems.jl DescriptorStateSpace object.

to_ode_problem(model::NonstationaryModel[, p=nothing; u=(t->0), x0=0.0, tspan=(0,1)])

Creates an ODEProblem for the model <: NonstationaryModel for a given input u(t). To use this method requires importing OrdinaryDiffEq as this functionality lives in an extension.

bode(model::LTIModel, ω::Real[, p=nothing; first=true])

Returns the transfer function evaluated at s=im*ω, C * (sE - A)^(-1) B + D, evaluated at the [1,1] entry if first==true.

bode(model::LTIModel, ωs::AbstractVector{<:Union{AbstractVector,Real}}[, p=nothing; first=true])

Returns the transfer function evaluated at s=im*ω for ω in ωs, evaluated at the [1,1] entry if first==true.

StabilityResidualErrorEstimator <: ErrorEstimator

Error estimator for a LinearModel which approximates the stability factor with stability_estimator(p) and computes the residual norm with residual_norm.

form_rom(pod_reductor, r=-1)

Pulls the first r left singular vectors from pod_reductor.decomp, and then Galerkin projects pod_reductor.model onto that basis, and returns the resulting ROM.

form_rom(sg_reductor, r=-1)

Calls galerkin_project on the FOM and returns a ROM with RB of dimension r. If r=-1, uses all available columns of sg_reductor.V.

form_rom(wg_reductor, r=-1)

Calls galerkin_project on the FOM and returns a ROM with RB of dimension r. If r=-1, uses all available columns of wg_reductor.V.

form_rom(bt_reductor[, r=-1])

Uses Petrov-Galerkin on the model to form a ROM of order r (largest possible if r==-1). Also, initializes it to bt_reductor.p if not nothing.

lift(pod_reductor, x_r)

Given a solution array x_r to a ROM formed by the pod_reductor lifts the solution(s) to the same dimension of the FOM.

lift(sg_reductor, x_r)

Given a solution array x_r to a ROM formed by the sg_reductor lifts the solution(s) to the same dimension of the FOM.

lift(wg_reductor, x_r)

Given a solution array x_r to a ROM formed by the wg_reductor lifts the solution(s) to the same dimension of the FOM.

lift(bt_reductor, x_r)

Given a solution array x_r to a ROM formed by the bt_reductor lifts the solution(s) to the same dimension of the FOM.

add_to_rb!(pod_reductor, snapshots[; noise=1])

Directly updates pod_reductor with new snapshots given in the columns of the matrix snapshots.

add_to_rb!(pod_reductor, parameters[; noise=0, progress=false])

Loops through the vector of parameters, forms their full order solutions, adds them to pod_reductor.snapshots, and then updates the singular values and singular vectors in pod_reductor.S and pod_reductor.V.

add_to_rb!(sg_reductor, snapshots[; noise=0])

Directly updates sg_reductor with new snapshots given in the columns of the matrix snapshots.

add_to_rb!(sg_reductor, parameters[; noise=0, progress=false])

Loops through the vector of parameters, forms their full order solutions, adds them to sg_reductor.snapshots, and then updates the reduced basis in sg_reductor.V.

add_to_rb!(wg_reductor, params[; noise=0, progress=false, eps=0.0, zero_tol=1e-15])

Loops through the vector of parameters params, computes the approximate estimator for each, selects the one with the highest error, and updates wg_reductor with the corresponding full order solution. Returns true if a vector is added to the RB, false otherwise.

add_to_rb!(wg_reductor, params, r[; noise=0, eps=0.0, zero_tol=1e-15])

Adds to wg_reductor at least r times by calling add_to_rb!(wg_reductor, params, noise=noise, eps=eps, zero_tol=zero_tol) several times. If all r are added, returns true, otherwise false.

get_rom(wg_reductor)

Helper method for getting the ROM from the wg_reductor object. Can otherwise obtain it through wg_reductor.rom.

AffineParametrizedArray(arrays, makeθi[, precompθ=false])

Struct for containing an affine parametrized array A(p) of the form A(p) = ∑ makeθi(p,i) * arrays[i].

Arrays are stored in the vector arrays for quick recomputation for new parameter values. Each must be of the same size and dimension so that they can be broadcast summed.

If precompθ set to false (default):

The function makeθi(p,i) takes in a parameter object first, and second an index i=1,...,length(arrays), and returns a scalar

If precompθ set to true:

The function makeθi(p) takes in a parameter object, and returns a vector of each of the affine terms.

Given aparr <: AffineParametrizedArray, and a new parameter value p, can form the full array A(p) by calling aparr(p).

formArray!(aparr, arr, p)

Given an array arr with the same dimensions as the arrays in the APArray aparr, form A(p) and place its values in arr.

eim(arrFun, param_disc[, ϵ=1e-2; maxM=100, noise=1])

Method for constructing an APArray object from a non-affinely parameter dependent matrix arrFun(p) by empirical interpolation.

param_disc must be a matrix with columns as parameter vectors, or a vector with elements as parameters.

Loops over the given parameter discretization until a maximum ∞-norm error of ϵ is achieved over the entire discretization, or until the maximum number of parameter values are chosen, given by maxM.

noise dictates the amount of printed output. Set 0 for no output, 1 for some output, ≥2 for most.

VectorOfVectors{T} <: AbstractMatrix{T}

Type for defining matrices as vectors of vectors for quick insertion of rows and columns. Stores the columns as vectors in vecs and its size in size.

Construct an empty VOV of dimensions (nrows × ncols) where one of nrows or ncols must be zero:

VectorOfVectors(nrows=0, ncols=0, T=Float64)

Construct a VOV from a Julia vector of vectors:

VectorOfVectors(vecs::AbstractVector{<:AbstractVector{T}})

Construct a VOV from a matrix:

VectorOfVectors(M::AbstractMatrix{T})

addRow!(vov::VectorOfVectors)

Add a row to the vector of vectors by appending zero(T) to each vector.

addRow!(vov::VectorOfVectors, row::AbstractVector)

Add the vector row to the last column of vov.

removeRow!(vov::VectorOfVectors)

Remove a row from the vector of vectors.

addCol!(vov::VectorOfVectors)

Add a column to the vector of vectors by appending zeros(T,nrow) to vov.vecs.

addCol!(vov::VectorOfVectors, col::AbstractVector)

Add the vector row to the last column of vov.

removeCol!(vov::VectorOfVectors)

Remove a column from the vector of vectors.

SCM{P} <: AbstractSCM <: Function where P <: Int

An SCM object for a SPD APArray matrix Ap. For non-SPD matrix Ap, can form Grammian matrix Ap'Ap and perform SCM on it.

scm = SCM(Ap::APArray, ps::AbstractVector[, ϵ=0.8, Mα=20, Mp=0; coercive=false, optimizer=HiGHS.Optimizer, lp_attrs=Dict(), make_monotonic=true, max_iter=500, noise=0])

Constructs an SCM object for Ap about parameter values ps. If coercive=false, then constructs the Grammian matrix Ap(p)'Ap(p). If hermitianpart(Ap(p)) is positive definite for all p, then set coercive=true. Constrains the SCM object to a relative gap of ϵ. If set ϵ=nothing, then does not perform constraining. Mα and Mp are the stability and positivity parameters, optimizer is the optimizer used for solving LPs in JuMP. lp_attrs can contain additional attributes to pass into the JuMP model (optimizer dependent). make_monotonic ensures that the lower-bound predictions increase monotonically. max_iter is the maximum number of SCM constraining iterations. noise determines the amount of printed output.

The scm object may be used by calling it on a parameter vector

scm(p[, Mα=20, Mp=0; which=:L, noise=0])

Which returns the following:

• which=:L - The lower-bound prediction at p
• which=:U - The upper-bound prediction at p
• which=:E - The relative gap 1 - LB/UB at p

ANLSCM{P} <: AbstractSCM <: Function where P <: Int

An SCM object for an APArray matrix Ap which uses nonlinear constraining.

scm = ANLSCM(Ap::APArray, ps::AbstractVector[, ϵ=0.8, Mα=20; optimizer=Ipopt.Optimizer, lp_attrs=Dict(), nonlin_alpha=0.0, max_iter=500, noise=0, constrain_kwargs...])

Constructs an ANLSCM object for Ap about parameter values ps. Constrains the ANLSCM object to a relative gap of ϵ. If set ϵ=nothing, then does not perform constraining. Mα is the stability parameter. optimizer is the optimizer used for solving NLPs in JuMP. lp_attrs can contain additional attributes to pass into the JuMP model (optimizer dependent). nonlin_alpha is the initial value of nonlinear constraint coefficients (0 for lightest constraint, 1 for tightest and LB guarantee). max_iter is the maximum number of SCM constraining iterations. noise determines the amount of printed output. Additional kwargs passed to constrain!.

The scm object may be used by calling it on a parameter vector

scm(p[, Mα=20; which=:L, noise=0])

Which returns the following:

• which=:L - The lower-bound prediction at p
• which=:U - The upper-bound prediction at p
• which=:E - The relative gap 1 - LB/UB at p

NNSCM{P} <: AbstractSCM <: Function where P <: Int

An SCM object for an APArray matrix Ap which uses domain decomposition and works in a natural norm. See Huynh et al., 2010, A natural norm...

scm = NNSCM(Ap::APArray, ps::AbstractVector[, ϵ=0.8, ϵ_β=0.8, ϕ=0.0, Mα=20; optimizer=HiGHS.Optimizer, lp_attrs=Dict(), max_iter=500, max_inner_iter=500, noise=0, constrain_kwargs...])

Constructs an NNSCM object for Ap about parameter values ps. Constrains the NNSCM object to a relative gap of ϵ. If set ϵ=nothing, then does not perform constraining. ϵ_β is the natural-norm relative gap ensured within subdomains. ϕ is the minimum natural-norm lower-bound prediction required to be considered within a domain. Mα is the stability parameter. optimizer is the optimizer used for solving NLPs in JuMP. lp_attrs can contain additional attributes to pass into the JuMP model (optimizer dependent). max_iter is the maximum number of SCM constraining iterations and max_inner_iter is the maximum number of iterations within a particular domain decomposition. noise determines the amount of printed output. Additional kwargs passed to constrain!.

The scm object may be used by calling it on a parameter vector

scm(p[, Mα=20; pbar=nothing, which=:L, noise=0])

Which returns the following:

• which=:L - The lower-bound prediction at p
• which=:U - The upper-bound prediction at p
• which=:E - The relative gap 1 - LB/UB at p

If pbar is not nothing, then if pbar in keys(scm.UBs), solves the above restricted to the domain decomposition about pbar.

copy_scm(scm::SCM[; lp_attrs=Dict(), noise=0])

Makes a copy of scm and provides lp_attrs to the new JuMP model.

copy_scm(scm::ANLSCM[; nonlin_alpha=0.0, lp_attrs=Dict(), noise=0])

Makes a copy of scm with initial nonlinear constraint coefficients nonlin_alpha and provides lp_attrs to the new JuMP model.

copy_scm(scm::NNSCM[; lp_attrs=Dict(), noise=0])

Makes a copy of scm and provides lp_attrs to the new JuMP model.

constrain!(scm::SCM, ps::AbstractVector, ϵ::Real, Mα::Int, Mp::Int[; make_monotonic=true, max_iter=500, noise=0, reigkwargs...])

Constrains scm about parameters ps to relative gap ϵ. Mα and Mp are the stability and positivity constraints respectively. make_monotonic ensures LB predictions increase monotonically. max_iter is the maximum number of SCM iterations ran. noise determines amount of printed output. See documentation of ModelOrderReductionToolkit.reig for arguments to reigkwargs.

constrain!(scm::ANLSCM, ps::AbstractVector, ϵ::Real, Mα::Int[; adaptive_nl_rate=1.1, adaptive_eps_tol=1e-6, max_iter=500, noise=0, reigkwargs...])

Constrains scm about parameters ps to relative gap ϵ. Mα is the stability constraint. adaptive_nl_rate and adaptive_eps_tol determine when to and rate at which to loosen nonlinear constraints. Set adaptive_nl_rate=nothing for no adaptive updating. max_iter is the maximum number of SCM iterations ran. noise determines amount of printed output. See ModelOrderReductionToolkit.reig for arguments to reigkwargs.

constrain!(scm::NNSCM, ps::AbstractVector, pbar::AbstractVector[, ϵ_β=0.8, ϕ=0.0, Mα=20, ps_left=trues(length(ps)); max_iter=500, p_choice=3, noise=0, reigkwargs...])

Constrains the pbar-domain decomposition of the scm about parameters ps. ϵ_β determines the natural-norm relative gap, ϕ determines the minimum LB prediction required for inclusion in domain, Mα is the stability parameter, ps_left is a vector of which parameters to consider, max_iter is maximum number of iterations.

p_choice can be 1, 2, or 3.

• p_choice=1 - Only add parameter with highest natural norm relative gap
• p_choice=2 - Only add parameter within domain (LB ≥ ϕ) with highest natural norm relative gap
• p_choice=3 - Do both of above

See ModelOrderReductionToolkit.reig for arguments to reigkwargs.

constrain!(scm::NNSCM, ps::AbstractVector, ϵ::Real[, ϵ_β=0.8, ϕ=0.0, Mα=20; pruning=false, removals=false, eps_keep=0.2, p_choice=3, max_iter=500, max_inner_iter=500, noise=0, reigkwargs...])

Constrains scm about parameters ps to relative gap ϵ with inner-maximum relative gap ϵ_β and domain-inclusion parameter ϕ. Mα is the stability constraint. pruning determines whether or not to exclude parameters from consideration when their relative gap decreases below ϵ for a greedier procedure. removals determines whether or not to try removing unnecessary domains at the start of each iteration with tolerance eps_keep. See other constrain! method for p_choice. max_iter is the maximum number of SCM iterations ran and max_inner_iter determines the maximum number of iterations per domain decomposition. noise determines amount of printed output.

Sigma_Min_RBF

A mutable struct used for approximating the minimum singular value of a matrix with parametric dependence A(p). Given a new parameter value p, approximate σ_min(A(p)) with sigma_min_rbf(p).

min_sigma_rbf(params, Ais, makeθAi[, ϕ=gaussian_rbf])

Method to form a interpolatory radial-basis function functor to approximate the minimum singular value of a parametrized matrix A(p). Pass as input either a matrix params where each column is a parameter vector, or a vector of vectors params where each vector is a parameter.

Optional argument for the radial-basis function ϕ, which defaults to the Gaussian ϕ(r) = exp(-r^2).

Returns a functor sigma_min_rbf <: Sigma_Min_RBF such that given a new parameter vector p, sigma_min_rbf(p) returns an approximation to the minimum singular value of A(p).

Offline, solves for the minimum singular value for each parameter in params, and uses this to form an interpolatory approximation in the form of log(σ_min(A(p))) ≈ ω_0 + ∑ (ω_i p_i) + ∑ (γ_i ϕ(p - p_i)) where p_i is the i'th parameter in params, and ω_0, ω_i, and γ_i are determined by the given params such that the approximation holds true for all p_i, and that ∑ γ_i = 0, and that ∑ γ_i p_i[j] = 0 for each j.

update_sigma_rbf!(sigma_min_rbf, ϕ)

Given a Sigma_Min_RBF object sigma_min_rbf, update its RBF, ϕ, and update the coefficients.

res <: ResidualNormComputer{T}

Abstract type for computing the norm of the residual for weak greedy RB methods. Must implement the update!(res, v) and compute(res, u_r, p) methods.

StandardResidualNormComputer{T} <: ResidualNormComputer{T} StandardResidualNormComputer(Ap::APArray,bp::APArray,V=nothing,X=nothing)

A struct for containing the necessary vectors and matrices for quickly compute the X-norm of the residual, r(u_r,p) = A(p) (u - V u_r) = b(p) - A(p) V u_r, by taking advantage of affine parameter dependence of A(p) and b(p).

Here, u solves A(p) u = b(p) with A and b having affine parameter dependence, and V is a matrix with columns defining bases for approximation spaces u ≈ V u_r.

ProjectionResidualNormComputer{T} <: ResidualNormComputer{T} ProjectionResidualNormComputer(Ap::APArray,bp::APArray,V=nothing,X=nothing)

A struct for containing the necessary vectors and matrices for quickly compute the X-norm of the residual, r(u_r,p) = A(p) (u - V u_r) = b(p) - A(p) V u_r, by taking advantage of affine parameter dependence of A(p) and b(p). Uses a projection method which is more stable than the standard method.

Here, u solves A(p) u = b(p) with A and b having affine parameter dependence, and V is a matrix with columns defining bases for approximation spaces u ≈ V u_r.

full_lu(A; steps=-1)

Performs a completely pivoted LU factorization on the matrix A, returning permutation vectors Q and P, and lower and upper triangular matrices L and U, such that if all steps are performed, then A[P,Q] = L*U.

reig(A::AbstractMatrix, [B=I; which=:L, kmaxiter=1000, noise=0, egwith=:arpack, restarts=100, krylovsteps=8, eps=1e-14, randdisks=1000, ignoreB=false, minstep=1.0, force_sigma=nothing])

Given (symmetric) matrices A and B, computes a real eigenvalue

$A x = λ B x$.

• which=:L corresponds to the largest eigenvalue
• which=:S corresponds to the smallest eigenvalue
• which=:SP corresponds to the smallest positive eigenvalue

Returns a tuple (λ,v) where λ<:Real is the eigenvalue and v the eigenvector.

Attempts to do this by shift-and-invert. Uses Arpack.jl if egwith=:arpack or ArnoldiMethod.jl if egwith=:arnoldimethod where the shifts are determined by the eigenvalue seeked and the Gershgorin disks of B⁻¹A.

If egwith=:arpack, kmaxiter determines maximum number of iterations. If egwith=:arnoldimethod, uses restarts to determine number of restarts.

NOTE: A randomized algorithm is used by default for determining Gershgorin disks. This is to speed up the case when size(A,1) ≫ 1.

Following paramers relate to selection of shifts:

• krylovsteps determines number of logarithmically spaced shifts are attempted
• ignoreB determines whether to approximate Gershgorin disks of B⁻¹A or A
• randdisks determines number of (random) columns of B⁻¹A or A to subselect

to compute Gershgorin disks of, others are ignored. Set randdisks ≥ size(A,1) for no randomization.

• minstep determines the minimum step size between logarithmically spaced shifts
• force_sigma forces trying shift and invert about this value if not nothing

Parameter eps is used to perturb shift parameters by relative amount to attempt to ensure no singular shifts. If the relative gap between the shift and the found eigenvalue is greater than reldifftol and the absolute gap is greater than absdifftol, attempts to re-solve by shifting about the found eigenvalue.

smallest_sval(A[; kmaxiter=1000, noise=0, reigkwargs...])

Given a matrix A, attempts to compute the smallest singular value of it by Krylov iteration and inversion. If unsuccessful, computes a full, dense svd. See docs for ModelOrderReductionToolkit.reig for reigkwargs options. Returns the smallest singular value, σ_min<:Real

largest_sval(A[; kmaxiter=1000, noise=0])

Given a matrix A, attempts to compute the largest singular value of it by Krylov iteration. If unsuccessful, computes a full, dense svd. Returns the largest singular value, σ_max<:Real.

orthonormalize_mgs2!(u, V)

Given a matrix V, and a new vector u, orthogonalize u with respect to the columns of V, and computes its norm nu. If nu != 0, divides u by nu and returns nu. If nu == 0, then u lives in the span of V, and 0 is returned.