Home Original page

emgr -- EMpirical GRamian Framework

Model Reduction Software

Name

emgr - EMpirical GRamian framework for (nonlinear) input-output systems.

Description

(image) emgr logo

System gramians are matrices associated to linear input-output systems and quantify their controllability, observability and minimality. Empirical system gramian matrices (empirical gramians) are a generalization computable for linear but also for parametric and nonlinear (state-space) systems and have, among others (see below), applications in model order reduction. Model reduction via empirical gramians is applicable to the state-space, to the parameter-space or both, through combined state and parameter reduction. For state reduction, the empirical reachability gramian / empirical controllability gramian and the empirical observability gramian, or the empirical cross gramian (the empirical linear cross gramian for large-scale linear systems) are available. For sensitivity analysis, parameter identification, and parameter reduction, the empirical sensitivity gramian (controllability of parameters), or the empirical identifiability gramian and empirical cross-identifiability gramian (observability of parameters) are provided. Combined reduction is enabled by the empirical joint gramian, which computes minimality of states (via the cross gramian) and observability of parameters (via the cross-identifiability gramian) simultaneously. Furthermore, the empirical output reachability gramian / empirical output controllability gramian is computable. Overall, empirical system Gramians are an advanced data-driven system-theoretic method and a numerical tool for GRAMIAN-based (non)linear input-output system analysis. The empirical gramian framework — emgr — is a compact open-source toolbox, compatible with OCTAVE and MATLAB, that provides a unified interface for unsupervised learning of system-theoretic operators from small data.

Scope

(image) emgr box
  • Model Order Reduction (MOR) | Model Reduction
    • Parametric Model Order Reduction (pMOR) | Robust Reduction
    • Nonlinear Model Order Reduction (nMOR)
    • Parameter Reduction
    • Combined State and Parameter Reduction (Combined Reduction)
  • Decentralized Control | Control Configuration Selection
  • Sensitivity Analysis
  • Parameter Identification | Input-Output Identifiability
  • Nonlinearity Quantification
  • Uncertainty Quantification
  • System Norms | System Indices | System Invariants
  • Optimal Sensor Placement | Optimal Actuator Placement
  • Matrix Equations

Download

Get emgr here: emgr.m (Version: 5.99, CRC32: f76f4638)    [emgr-5.99.zip]    [mirror]    [source]    [meta]    [repo]

(emgr is written in the Matlab programming language and requires Octave >= 5.2 or MATLAB >= 2017b, yet emgr has no dependencies on toolboxes or packages.)

License

All source code is licensed under the open source BSD-2-Clause license:

Copyright (c) 2013–2022, Christian Himpe

All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Disclaimer

This is research software!

Model

The target mathematical model is a general input-output system (aka control system): ẋ(t) = f(x(t),u(t),p,t) y(t) = g(x(t),u(t),p,t) consisting of a first-order ordinary differential equation (ODE) with vector field f and an algebraic output function g, of which either (may) depend on input u(t), state x(t), parameter p and time t, yielding the output y(t). This model maps an input (control) to an output (quantities of interest) via the state evolving under a dynamical system.

Algorithm

The computation of empirical system Gramians is described in the open-access article:

   C. Himpe; "emgr -- the Empirical Gramian Framework"; Algorithms, 11(7): 91, 2018;

and references therein.

Usage

General Usage: W = emgr(f,g,s,t,w,pr,nf,ut,us,xs,um,xm,dp);

Minimal Usage: W = emgr(f,g,s,t,w);

About Info Usage: v = emgr('version');

Mandatory Arguments

  • f - handle to a function with signature x = f(x,u,p,t) , the system's vector field,
  • g - handle to a function with signature y = g(x,u,p,t) , the system's output functional; g = 1 implies y = x,
  • s - three component vector s = [M,N,Q] holding number of inputs (M), states (N) and outputs (Q),
  • t - two component vector t = [dt,Tf] holding time step (dt) width and time horizon (Tf),
  • w - a character selecting the gramian type (For details see Gramians).

Optional Arguments

(flowchart) gramian interrelations
  • pr - System's parameters (Default: 0; 's', 'i', 'j' require two columns for min and max parameter),
    • vector - column vector holding the parameters,
    • matrix - set of column vectors holding different sets of parameters.
  • nf - Thirteen component vector holding options; (Default: 0, for details see Option Flags).
  • ut - Input time series (Default: 'i'),
    • handle - handle to function with signature u_t = ut(t),
    • 'i' - delta impulse input,
    • 's' - step input / load vector / source term,
    • 'h' - decaying exponential chirp input via havercosine,
    • 'a' - cardinal sine / sinc input,
    • 'r' - uniformly random binary input.
  • us - Steady-state input (Default: 0),
    • scalar - sets all M input components to provided value,
    • vector - steady-state input column vector of dimension M.
  • xs - Steady-state, also nominal initial state (Default: 0),
    • scalar - sets all N steady-state components to provided value,
    • vector - steady-state column vector of dimension N.
  • um - Input perturbation scales (Default: 1),
    • scalar - set all maximum scales to argument,
    • vector - set maximum to scale to argument,
    • matrix - set scales to argument, used as is.
  • xm - Initial state perturbation scales (Default: 1),
    • scalar - set all maximum scales to argument,
    • vector - set maximum to scale to argument,
    • matrix - set scales to argument, used as is.
  • dp - Custom inner product handle xy = dp(x,y) (Default: @mtimes).

Gramian Types

(chart) #emgr application matrix
  • 'c' - Empirical Controllability Gramian (WC),
  • 'o' - Empirical Observability Gramian (WO),
  • 'x' - Empirical Cross Gramian (WX sometimes noted as WCO),
  • 'y' - Empirical Linear Cross Gramian (WY),
  • 's' - Empirical Sensitivity Gramian (WS, controllability of parameters),
  • 'i' - Empirical Identifiability Gramian (WI, observability of parameters, via augmented observability Gramian),
  • 'j' - Empirical Joint Gramian (WJ, observability of parameters via cross-identifiability Gramian).

Option Flags

  • nf(1) - Center time series around (Default: 0),
    • = 0 - zero (balanced POD),
    • = 1 - steady-state (empirical covariances),
    • = 2 - final state,
    • = 3 - arithmetic average over time (empirical gramians),
    • = 4 - root-mean-square over time,
    • = 5 - mid-range over time.
  • nf(2) - Input scale sequence (Default: 0),
    • = 0 - single um = um;
    • = 1 - linear um = um * [0.25, 0.50, 0.75, 1.0];
    • = 2 - geometric um = um * [0.125, 0.25, 0.5, 1.0];
    • = 3 - logarithmic um = um * [0.001, 0.01, 0.1, 1.0];
    • = 4 - sparse um = um * [0.01, 0.5, 0.99, 1.0];
  • nf(3) - State scale sequence (Default: 0),
    • = 0 - single xm = xm;
    • = 1 - linear xm = xm * [0.25, 0.50, 0.75, 1.0];
    • = 2 - geometric xm = xm * [0.125, 0.25, 0.5, 1.0];
    • = 3 - logarithmic xm = xm * [0.001, 0.01, 0.1, 1.0];
    • = 4 - sparse xm = xm * [0.01, 0.5, 0.99, 1.0];
  • nf(4) - Input transformations (Default: 0),
    • = 0 - unit um = [-um, um];
    • = 1 - single um = um;
  • nf(5) - State transformations (Default: 0),
    • = 0 - unit xm = [-xm, xm];
    • = 1 - single xm = xm;
  • nf(6) - Gramian normalization, only WC, WO, WX, WY (Default: 0),
    • = 0 - no normalization,
    • = 1 - normalize with Jacobi preconditioner,
    • = 2 - normalize with steady-state (input).
  • nf(7) - State gramian variant (Default: 0),
    • = 0 - regular,
    • = 1 - output-controllability gramian (WC, WS),
    • = 1 - averaged observability gramian (WO, WI),
    • = 1 - non-symmetric cross gramian (cross operator WZ) for non-square, non-symmetric or non-gradient systems (WX, WY, WJ).
  • nf(8) - Enable nominal input during parameter and state perturbations, only WO, WX, WS, WI, WJ (Default: 0),
    • = 0 - no extra input,
    • = 1 - extra input for state or parameter perturbations.
  • nf(9) - Parameter centering and scale sequence, only WS, WI, WJ (Default: 0),
    • = 0 - no centering and linear scales,
    • = 1 - linear centering and scales,
    • = 2 - logarithmic centering and scales.
    • = 3 - nominal centering and scales.
  • nf(10) - Parameter gramian variant, only WS, WI, WJ (Default: 0),
    • = 0 - input-state average (WS),
    • = 1 - input-output average (WS),
    • = 0 - coarse Schur-complement (WI, WJ),
    • = 1 - approximate Schur-complement (WI, WJ).
    • = 2 - exact Schur-complement (WI, WJ).
  • nf(11) - Cross gramian partition width, only WX, WJ (Default: 0),
    • = 0 - full cross Gramian, no partitioning,
    • < N - maximum partition width.
  • nf(12) - Partitioned cross gramian running index, only WX, WJ (Default: 0),
    • = 0 - no partitioning,
    • > 0 - index of cross gramian partition to be computed.
  • nf(13) - Weight time series (Default: 0),
    • = 0 - one (no weighting),
    • = 1 - linear time-weighting,
    • = 2 - quadratic time-weighting,
    • = 3 - state-based weighting,
    • = 4 - scale-based weighting.
    • = 5 - reciprocal square-root time-weighting.

Return Values

  • State-space empirical gramian matrix (for: WC, WO, WX, WY).
  • Cell array of state-space and parameter-space empirical gramian matrix (for: WS, WI, WJ).

Solver Interface

SSP32 Runge-Kutta schematic A custom solver for quantities of interest of an ordinary differential equations is set by a function handle in the global variable: ODE.
  • Function signature: y = solver(f,g,t,x0,u,p)
    • f - handle to a function with signature x = f(x,u,p,t) , the system's vector field,
    • g - handle to a function with signature y = g(x,u,p,t) , the system's output functional,
    • t - two component vector t = [dt,Tf] holding time step width (dt) and time horizon (Tf),
    • x0 - column vector of dimension N holding initial state,
    • u - handle to function with signature u_t = u(t),
    • p - column vector holding (current) parameter(s).
  • Included default solver: SSPx2

Minimal Example

A = -0.5*eye(4) % System matrix B = [0;1;0;1] % Input matrix C = [0,0,1,1] % Output matrix WX = emgr(@(x,u,p,t) A*x + B*u, @(x,u,p,t) C*x,[1,4,1],[0.1,10.0],'x') % ≈ B*C

Tests

Tests are defined in emgrTest.m and are evaluated using est. The tests are performed on a time-invariant, linear, state-space symmetric MIMO system and optional linear parametrization: ẋ = A*x + B*u + F*p y = C*x

Demos

The demo codes can be found in estDemo.m and are evaluated using est.

Combined Reduction: Nonlinear System (WJ)

Here a parametrized nonlinear control system is reduced in states and parameters using the empirical joint gramian. The nonlinearity is given by a hyperbolic tangent with parametrized argument. This is a parametrized hyperbolic network model from "Modeling and control based on a new neural network model" by Y. Quan et al (2001). ẋ = A*tanh(p.*x) + B*u y = C*x

(plot) combined_wj reduction #emgr
Relative L2 output error for varying state and parameter dimension.

Benchmark: Inverse Sylvester Procedure (WC+WO)

The inverse Sylvester procedure can be used to generate random linear systems based on "On generating random systems: a gramian approach" by S.C. Smith and J. Fisher (2003). Here, the required solution of the Sylvester equation is accomplished by the empirical gramians after factoring (-1). The generated random system is then reduced by using balanced truncation. 0 = A*WC + WC*A' - B'*B = (-1*a)*WC + WC*(-1*a') - B'*B = a*(-1*WC) + (-1*WC)*a' - B'*B

(plot) benchmark_ilp reduction #emgr
Relative L1, L2, L∞, L0 output error for varying state dimension.

Benchmark: Flexible Space Structures (WX)

A procedural benchmark modelling modal structural dynamics for applications such as truss masts for space environments (i.e.: COFS-1) based on "Model Reduction for Flexible Space Structures" by W. Gawronski and T. Williams (1991). This is a linear second order model, which is reduced by approximate balancing via the empirical cross gramian.

(plot) benchmark_fss #emgr
Relative L1, L2, L∞, L0 output error for varying state dimension.

Benchmark: Nonlinear Model Reduction (WX)

An implementation of a model reduction benchmark from "Model Reduction for Nonlinear Systems" by Y. Chen (1999). Here the nonlinear RC-ladder benchmark, a resistor-capacitor cascade with nonlinear resistors (diodes), is tested. The empirical cross gramian is utilized to reduce this nonlinear system.

(plot) benchmark_non #emgr
Relative L1, L2, L∞, L0 output error for varying state dimension.

Quadratic Output: Kernel Matrix (WX)

The randomly generated example from "Balanced Model Reduction via the Proper Orthogonal Decomposition" by K. Willcox and J. Peraire (2002), is used with a quadratic output functional, and reduced by a empirical cross kernel matrix utilizing a Gaussian kernel. ẋ = A*x + B*u y = x'*x

(plot) reduced quadratic output #emgr
Relative L1, L2, L∞, L0 output error for varying state dimension.

DMD-Galerkin: All-Pass System (WC)

An all pass system based on "Asymptotically Stable All-Pass Transfer Functions: Canonical Form, Parametrization and Realization" by R.J. Ober (1987), which has a single Hankel singular value of multiplicity of the system's order, is reduced in state dimension by the DMD-Galerkin method via the empirical controllability Gramian! ẋ = A*x + B*u y = C*x

(plot) reduced all pass #emgr
Relative L1, L2, L∞, L0 output error for varying state dimension.

PDE Reduction: Advection Equation (WY)

Discretized with a first order upwind finite difference scheme, the 1D linear transport equation is a hyperbolic partial differential equation, whereas the left boundary of the domain is set as input and the right boundary is set as output; as presented in "Efficient Calculation of Thermoacoustic Modes Utilizing State-Space Models" by M. Meindl, T. Emmert and W. Polifke (2016) as flame state-space model. This system is reduced by the empirical linear cross gramian and empirical dominant subspaces. dz/dt = -a dz/dx z(0,t) = u(t) y(t) = z(1,t)

(plot) reduced advection #emgr
Relative L1, L2, L∞, L0 output error for varying state dimension.

Nonlinear Second Order Reduction: 5-body Choreography (WO)

A 5-body choreography, based on the gravitational N-body problem, resembling the figure eight is reduced in states; based on "New Families of Solutions in N-Body Problems" by C. Simo (2001). This is a particularly difficult model to reduce.

(plot) reduced nbody #emgr
Reduced order 5-body figure-8 orbits.

Sensitivity Analysis: Stable Orbits Inside Black Holes (WS)

This is an implementation of: "Is there life inside black holes?" by V. Dokuchaev (2011), describing stable orbits of photons and particles inside the event horizon of a black hole. The underlying system is nonlinear and the output function is a coordinate transformation from Boyer-Lindquist to Cartesian coordinates. In the computation of particle and photon orbits the associated parameters are identified by the empirical sensitivity gramian.

(plot) blackhole orbit #emgr
Quasi-stable orbits inside the event horizon of a black hole.

About

A gram matrix or gramian matrix W is the result of all inner products of a set of vectors V = [v1 ... vn], in other words: W = VT V, and closely related to covariance matrix as well as cross-covariance matrix, correlation matrix or cross-correlation matrix respectively. Properties of (linear) control systems can be assessed by the system gramians, which are based on the controllability and observability operators. Classically, the controllability gramian and observability gramian are utilized in balancing model reduction methods. The cross gramian is not gramian matrix, but encodes controllability and observability information, and thus minimality, into a single matrix. Overall, system gramians quantifiy the system-theoretic properties controllability (reachability), observability and (structural) identifiability.

Empirical gramians extend this approach to nonlinear control systems and thus enable nonlinear model reduction. For linear systems, these empirical system gramians are equal to the classic system gramians, yet, empirical gramians are also computable for parametric and nonlinear systems. Furthermore, empirical gramians contain more information about the underlying system; and particularly the empirical cross gramian conveys even additional information. This makes empirical system gramians a versatile tool for mathematical system theory, control theory or computational science and engineering (CSE).

The (discrete) empirical cross gramian encloses information on the input-output behavior of the associated control system as well as approximate Hankel Singular Values and can be computed very efficiently. For large-scale linear systems, the linear empirical cross gramian, related to balanced POD, can be utilized. And for parametric systems, the (empirical) joint gramian, derived from the cross gramian, is available for combined reduction. In case of custom training input, instead of an empirical gramian matrix, an empirical covariance matrix can also be computed.

Model order reduction or short model reduction is an active research field in scientific computing and CSE, which investigates the algorithmic construction of low-order surrogate models for high-dimensional differential equation models. The resulting reduced order models allow a substantially more economical evaluation in terms of computational or memory resources compared to the original large-scale full order model. This numerical model reduction is essential, as developments in computing technology are constantly outpaced by the requirements for model simulation in dimensionality or complexity.

A primary application of empirical Gramians is data-driven model reduction of nonlinear input-output systems. In the language of data science, a nonlinear input-output system corresponds to a recurrent neural network, i.e. from (data-driven) system identification, which can also be (empirical) Gramian-based, or mathematical modeling. Empirical-Gramian-based model reduction then translates to unsupervised learning, where the empirical Gramians encode the network's reducibility based on dynamic behavior induced by systematic perturbations, correlating to synthetic data, and based on a model making it physics-informed.

References

Contact

Send feedback to: ch@gramian.de

Cite

gramian.de quick response code
  • Cite as: C. Himpe (2022). emgr — EMpirical GRamian Framework (Version 5.99) [Software]. https://gramian.de doi:10.5281/zenodo.6457616
  • BibTeX: @MISC{emgr, author={C.~Himpe}, title={{emgr -- EMpirical GRamian Framework} (Version~5.99)}, howpublished={\url{https://gramian.de}}, year={2022}, doi={10.5281/zenodo.6457616}}
  • DOI: 10.5281/zenodo.6457616 (Version 5.99)
  • Except where otherwise noted, content on this site is licensed under a CC BY 4.0 license.
  • Last Change: 2023-09-20

Links

Notes

  • emgr has AUTHORS, README, CHANGELOG, CITATION.cff, CODE, LICENSE, VERSION, RUNME.m files.
  • emgr has a special version emgr.py for Python 3.X.
  • emgr can compute the empirical cross gramian columnwise, and thus in parallel on distributed memory systems.
  • emgr is not explicitly parallelized but multi-core ready by extensive vectorization and implicit parallelization.
  • emgr has highlighted loops that qualify for explicit parallelization using parfor.
  • emgr's custom dot-product can also be used for GPGPU based matrix multiplication.
  • emgr core function consists of a single file and has less than 600 lines of code!
  • emgr's empirical gramian quality depends heavily on the specification of the operating region and the solver.
  • emgr's default solver can be used externally by calling emgr('version'), which sets the global variable ODE to a handle.

Troubleshooting

  • Issue: An empirical gramian contains very large or infinity values.
    • Fix: Likely an unstable trajectory, ensure the perturbations (scales) or parameters do not destabilize the system and a suitable ODE solver is used.
  • Issue: The linear cross gramian WY is zero.
    • Fix: The second argument for the empirical inear cross gramian has to be the system's adjoint vector field, NOT the output functional.
  • Issue: The identifiability gramian WI{2} or cross-identifiability gramian WJ{2} are zero.
    • Fix: The parameter steady state (initial state) does not excite the system; usually this means setting xs ≠ 0 or nf(8) = 1.
  • Issue: Using parfor results in a loop-variable error.
    • Fix: Since parfor cannot handle non-consecutive loop indices, zero scales or zero min/max parameters are not admissible. Constant parameters may be used if they are last in the parameter vector.
  • Issue: Using nf(10) = 1 yields infinity values in a parameter gramian ('s', 'i', 'j').
    • Fix: The min or max parameter sets likely contain zeros; try shifting parameters.
  • Issue: Using nf(1) = 1 with gramians 'o', 'x', 'i', 'j' results in wrong results.
    • Fix: In these configurations the steady-state output ys = g(xs,us,p(:,k),0) (for all k) is computed, check if g works accordingly.
  • Issue: The joint gramian / cross-identifiability gramian 'j' results in an error when using a Gaussian kernel as custom inner product.
    • Fix: Due to the memory-efficient computation via non-square intermediate matrices this combination (as well as with any radial kernel) is incompatible.

Frequently Asked Questions

  • Q: What is model order reduction? What is model reduction?
    • A: Model order reduction, or short model reduction, is the algorithmic computation of accurate low-dimensional surrogate models from high-dimensional differential-equation models.
  • Q: What are system Gramians?
    • A: System gramians are matrices characterizing the associated (linear) input-output system.
  • Q: What are empirical Gramians?
    • A: Empirical system Gramians, or short empirical Gramians, are data-driven approximate system Gramians for parametric or nonlinear systems.
  • Q: How are empirical Gramians related to machine learning?
    • A: In the language of machine learning, empirical Gramians provide physics-informed unsupervised learning with synthetic data on recurrent neural networks.
  • Q: What is emgr?
    • A: The empirical Gramian framework is a backend toolbox to compute empirical system Gramians.
  • Q: Why not just linearize a nonlinear input-output system?
    • A: Linearization typically uses a single linearization point, while empirical Gramians practically average multiple linearization points in a custom operating region of the system.
  • Q: Which files do I need for my own project from the emgr.zip archive?
    • A: Only the emgr.m file is needed; all other files are sample code, tests or documentation.

See Also

  • Model Reduction Routines (Another Empirical Gramian Software; empirical WC and WO only)
  • gram (MATLAB Control System Toolbox; linear WC and WO only)
  • gram (Octave Control Package; linear WC and WO only)

Category

(plot) cross gramian surface #emgr
  • MSC2020: 93-04
  • MSC2010: 93A15, 93B11, 93C10
  • ACM1998: G.4
  • Controllability, Observability
  • Science/Math/Software/MATLAB
  • Programming Languages: Matlab
  • math.OC, cs.SY, cs.MS
  • Topic :: Scientific/Engineering :: Mathematics
  • MSC2090: 285S29 (😛)

tl;dr

emgr: Unsupervised learning of I/O system properties for data-driven control.