## Model Reduction Software

### Name

emgr - **EM**pirical **GR**amian framework for (nonlinear) input-output systems.

### Description

System gramians are matrices associated to linear input-output systems and quantify their controllability, observability and minimality.
Empirical system gramian matrices (**empirical gramians**) are computable for linear but also nonlinear (state-space) systems and have, among others, application in model order reduction (MOR) or uncertainty quantification (UQ).
Model reduction using empirical gramians can be applied to the state-space, to the parameter-space or both through combined reduction.
For state reduction, empirical balanced truncation by the **empirical controllability gramian** and the **empirical observability gramian**, or alternatively, direct truncation (approximate balancing) by the **empirical cross gramian** (or the **empirical linear cross gramian** for large-scale linear systems) is available.
For parameter reduction, parameter identification and sensitivity analysis by the **empirical sensitivity gramian** (controllability of parameters) or the **empirical identifiability gramian** (observability of parameters) are provided.
Combined state and parameter reduction is enabled by the **empirical joint gramian**, which computes minimality of states (cross gramian) and observability of parameters (**empirical cross-identifiability gramian**) concurrently.
The empirical gramian framework - `emgr`

- is a compact open-source toolbox for GRAMIAN-based model reduction and compatible with OCTAVE and MATLAB.
This mathematical software provides a common interface for the computation of empirical gramians and empirical covariance matrices.

### Scope

- 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
- Sensitivity Analysis
- Structural Identifiability | Parameter Identification
- Optimal Sensor Placement | Optimal Actuator Placement
- Nonlinearity Quantification
- System Indices | System Invariants

- for Input-Output Systems | Control Systems | Dynamical Systems:
- linear & nonlinear,
- time-invariant & time-varying,
- dense & sparse,
- parametrized | parametric,
- ODE | spatially discretized PDE,
- control-affine

- Mathematical system (time t, state x, input u, parameter p)
- Vector Field f: ẋ(t) = f(x(t),u(t),p,t)
- Output Functional g: y(t) = g(x(t),u(t),p,t)

### Download

####
Get `emgr`

here: emgr.m (Version: 5.6, CRC32: 999979f4)
[mirror]
[source]
[meta]
[emgr_oct.m]
[emgr_lgc.m]
[emgr.py]

^{(emgr is written in the Matlab programming language and requires Octave or MATLAB, but 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--2019, 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!**

### 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 output functional;`g = 1`

implies`y = x`

,`s`

- three component vector`s = [M,N,Q]`

holding number of inputs, states and outputs,`t`

- two component vector`t = [h,T]`

holding time step width and time horizon,`w`

- a character selecting the gramian type (For details see Gramians).

### Optional Arguments

`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`

- Twelve component vector holding options; (Default:`0`

, for details see Option Flags).`ut`

- Input time series (Default:`1`

),`handle`

- handle to function with signature`u_t = ut(t)`

,`0`

- pseudo-random binary input,`1`

- delta impulse input,`∞`

- decaying exponential chirp input via havercosine.

`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 used as 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`z = dp(x,y)`

(Default:`@mtimes`

).

### Empirical Gramian Types

`'c'`

- Empirical Controllability Gramian (`W`

),_{C}`'o'`

- Empirical Observability Gramian (`W`

),_{O}`'x'`

- Empirical Cross Gramian (`W`

sometimes also noted as_{X}`W`

or_{CO}`X`

),_{CG}`'y'`

- Empirical Linear Cross Gramian (`W`

),_{Y}`'s'`

- Empirical Sensitivity Gramian (`W`

),_{S}`'i'`

- Empirical Identifiability Gramian (`W`

via augmented observability Gramian),_{I}`'j'`

- Empirical Joint Gramian (`W`

has cross-identifiability Gramian)._{J}

`W`_{Z}

can be selected by setting option flag `nf(7) = 1`

and is compatible with the empirical cross gramian `W`_{X}

, the empirical linear cross gramian `W`_{Y}

and the empirical joint gramian `W`_{J}

.
### Option Flags

`nf(1)`

- Center timer 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`

- midrange 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)`

- Normalizing (Default:`0`

),`= 0`

- no normalization,`= 1`

- normalize with Jacobi preconditioner,`= 2`

- normalize with steady-state (input).

`nf(7)`

- Non-Symmetric cross gramian, only`W`

(Default:_{X}, W_{Y}, W_{J}`0`

),`= 0`

- regular cross gramian,`= 1`

- non-symmetric cross gramian (cross operator`W`

) for non-square, non-symmetric or non-gradient systems._{Z}

`nf(8)`

- Enable nominal input during parameter and state perturbations, only`W`

(Default:_{O}, W_{X}, W_{S}, W_{I}, W_{J}`0`

),`= 0`

- no extra input,`= 1`

- extra input for state or parameter perturbations.

`nf(9)`

- Parameter centering and scale sequence, only`W`

(Default:_{S}, W_{I}, W_{J}`0`

),`= 0`

- no centering and linear scales,`= 1`

- linear centering and scales,`= 2`

- logarithmic centering and scales.

`nf(10)`

- Parameter gramian variant, only`W`

(Default:_{S}, W_{I}, W_{J}`0`

),`= 0`

- input-state average (`W`

), detailed Schur-complement (_{S}`W`

),_{I}, W_{J}`= 1`

- input-output average (`W`

), approximate Schur-complement (_{S}`W`

)._{I}, W_{J}

`nf(11)`

- Cross gramian partition width, only`W`

(Default:_{X}, W_{J}`0`

),`= 0`

- full cross Gramian, no partitioning,`< N`

- maximum partition width.

`nf(12)`

- Partitioned cross gramian running index, only`W`

(Default:_{X}, W_{J}`0`

),`= 0`

- no partitioning,`> 0`

- index of cross gramian partition to be computed.

### Return Values

- State-space empirical gramian matrix (for:
`W`

)._{C}, W_{O}, W_{X}, W_{Y} - Cell array of state-space and parameter-space empirical gramian matrix (for:
`W`

)._{S}, W_{I}, W_{J}

### Configuration

Custom solver (via handle in global variable:`ODE`

) for quantity of interest of an ordinary differential equations.
- 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 output functional,`t`

- two component vector`t = [h,T]`

holding time step width and time horizon,`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
- Explicit Second Order Runge-Kutta (RK2)
- Optimal Strong Stability Preserving (SSP)
- Configurable number of stages for improved stability (Default: 3)

### Utility

The platform framework curios.m (Clearing Up Reducibility of Input-Output Systems) enables testing and evaluating methods based on empirical gramians.`curios`

has the following capabilities:
- Model Reduction (controllability-truncation, observability-truncation, linear-direct-truncation, nonlinear-direct-truncation, linear-balanced-truncation, nonlinear-balanced-truncation, linear-dominant-subspaces, nonlinear-dominant-subspaces),
- Parameter Reduction (controllability-based, observability-based, minimality-based),
- Combined Reduction (controllability-based, observability-based, minimality-based),
- Sensitivity Analysis (input-state-based, input-output-based),
- Parameter Identification (state-output-based, input-output-based),
- Decentralized Control (linear, nonlinear),
- Nonlinearity Quantification (input-based, state-based, output-based),
- State Index (controllability, observability, minimality),
- System Index: (cauchy-index, system-entropy, system-gain, hinf-bound, hankel-bound, system-symmetry, energy-fraction, storage-efficiency, nyquist-area, robustness-index, recoverability-index, input-output-coherence).

### Minimal Code 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`curios`

.
The tests are performed on a time-invariant, linear, state-space symmetric MIMO system, with a negative Lehmer matrix `A`

as a system matrix and optional linear parametrization:
```
ẋ = A*x + B*u + F*p
y = C*x
```

### Demos

The demo codes can be found in examples.m and are evaluated using`curios`

.
## Combined Reduction: Nonlinear System (
_{J} | |

## Benchmark: Inverse Sylvester Procedure (
_{C}+W_{O} | |

## Benchmark: Flexible Space Structures (
_{X} | |

## Benchmark: Nonlinear Model Reduction (
_{X} | |

## PDE Reduction: Advection Equation (
_{Y} | |

## Nonlinear Second Order Reduction: 5-body Choreography (
_{O} | |

## Sensitivity Analysis: Stable Orbits Inside Black Holes (
_{S} |

### About

A **gramian matrix** `W`

is the result of all inner products of a set of vectors `V = [v1 ... vn]`

, in other words: `W = V`

.
Properties of (linear) control systems can be assessed by the ^{T} V**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 encodes controllability and observability information into a single matrix and thus minimality.
Overall, empirical 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 gramians are equal to the classic system gramians,
yet empirical gramians also computable for parametric and nonlinear systems.
Furthermore, empirical gramians contain more information about the underlying system;
and the empirical cross gramian conveys even additional information.
This makes empirical system gramians a versatile tool for mathematical engineering.

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 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 input, an empirical covariance matrix can also be computed.

**Model reduction** or **model order reduction** is an active research field in *scientific computing* and *computational science and engineering*,
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.
Model reduction is essential as developments in computing technology are constantly outpaced by the requirements for model simulation in dimensionality or complexity.

### References

- P. Benner, C. Himpe; "Cross-Gramian-Based Dominant Subspaces"; arXiv, math.OC: 1809.08066, 2018.
- C. Himpe; "emgr - the Empirical Gramian Framework"; Algorithms, 11(7): 91, 2018.
*Synopsis:*Mathematical and technical documentation of the empirical Gramian framework detailing all features and capabilities.*Links:*self-archived, source, runmycode

- C. Himpe, T. Leibner, S. Rave, J. Saak; "Fast Low-Rank Empirical Cross Gramians"; in: Proceedings in Applied Mathematics and Mechanics, 17: 841--842, 2017.
- C. Himpe, M. Ohlberger; "Cross-Gramian-Based Model Reduction: A Comparison"; in: Model Reduction of Parametrized Systems: 271--283, 2017.
*Synopsis:*Summary and tests for Sylvester-equation-based cross gramian, empirical linear cross gramian and empirical cross gramian.*Links:*preprint, source-code, runmycode

- C. Himpe; "Combined State and Parameter Reduction for Nonlinear Systems with an Application in Neuroscience"; Westfälische Wilhelms Universität, Sierke Verlag Göttingen, 2017.
*Synopsis:*Combined state and parameter reduction via empirical Gramians.*Links:*self-archived, source-code, zenodo

- C. Himpe, M. Ohlberger; "A note on the cross gramian for non-symmetric systems"; System Science and Control Engineering, 4(1): 199--208, 2016.
*Synopsis:*A new cross gramian matrix for non-symmetric, non-gradient and non-square system.*Links:*self-archived, source-code, runmycode

- C. Himpe, M. Ohlberger; "Accelerating the Computation of Empirical Gramians and Related Methods"; at: 5th International Workshop on Model Reduction in Reacting Flows, 2015.
*Synopsis:*Numerical performance enhancements for empirical gramians.*Links:*self-archived, source-code

- C. Himpe, M. Ohlberger; "The Empirical Cross Gramian for Parametrized Nonlinear Systems"; in: IFAC-PapersOnLine (Vienna International Conference on Mathematical Modelling), 48(1): 727--728, 2015.
*Synopsis:*Use of the empirical cross gramian for parametric model order reduction by averaging.*Links:*source-code, runmycode

- C. Himpe, M. Ohlberger; "Cross-Gramian-Based Combined State and Parameter Reduction for Large-Scale Control Systems"; Mathematical Problems in Engineering, 2014: 843869, 2014.
*Synopsis:*Introduction of the empirical cross gramian and the derived joint gramian as well as gramian-based combined reduction.*Links:*self-archived, source-code, runmycode

- C. Himpe, M. Ohlberger; "Model Reduction for Complex Hyperbolic Networks"; in: Proceedings of the European Control Conference: 2739--2743, 2014.
*Synopsis:*Cross gramian based model reduction of a time-varying non-symmetric system with application in network cosmology.*Links:*preprint, source-code, runmycode

- C. Himpe, M. Ohlberger; "A Unified Software Framework for Empirical Gramians"; Journal of Mathematics, 2013: 365909, 2013.
*Synopsis:*Mathematical background on and implementation of empirical gramians under a common interface.*Links:*self-archived, source-code, runmycode

### Contact

Send feedback to: `ch@gramian.de`

### Cite

- Cite as: C. Himpe (2019). emgr - EMpirical GRamian Framework (Version 5.6) [Software].
**https://gramian.de****doi:10.5281/zenodo.2530021** - BibTeX:
`@MISC{emgr, author={C.~Himpe}, title={{emgr - EMpirical GRamian Framework} (Version~5.6)}, howpublished={\url{https://gramian.de}}, year={2019}, doi={10.5281/zenodo.2530021}}`

- DOI: 10.5281/zenodo.2530021 (Version 5.6)
- Except where otherwise noted, content on this site is licensed under a CC BY 4.0 license.
- Last Change:
*2019-01-02*

### Links

- emgr-ref.pdf (emgr reference card)
- emgr-flyer.pdf (emgr flyer)
- emgr-poster-2018.pdf (emgr poster MATLAB Expo)
- emgr-poster-2015.pdf (emgr poster OctConf)
- github.com/gramian/emgr (emgr git repository at github)
- orms.mfo.de/project?id=345 (emgr at Oberwolfach References on Mathematical Software)
- researchgate.net/project/emgr-EMpirical-GRamian-Framework (emgr at researchgate)
- mathworks.com/matlabcentral/fileexchange/40169-emgr-empirical-gramian-framework (emgr at fileexchange)
- swmath.org/software/7554 (emgr at swMATH)
- freshcode.club/projects/emgr (emgr at freshcode)
- openhub.net/p/emgr (emgr at openhub)
- modelreduction.org (Model Order Reduction Wiki)
- morepas.org (Model Reduction for Parametrized Systems)
- git.io/hapod (Hierarchical Approximate Proper Orthogonal Decomposition)
- git.io/mtips (Octave / Matlab Snippets)

### Notes

**emgr**has a README.**emgr**has a special version emgr_oct.m using Octave's advanced syntax.**emgr**has a special version emgr_lgc.m for MATLAB 2016a and earlier.**emgr**has an experimental version emgr.py for Python 2.X and 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 be used for GPGPU based matrix multiplication.**emgr**consists of a single file and has less than 500 lines of code!**emgr**can only handle real-valued trajectory data.**emgr**'s empirical gramian quality depends heavily on the specification of the operating region and the solver.**emgr**can use low-dimensional steady-state perturbation by reconstructing the high-dimensional initial state via an encapsulated solver.

### 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`W`

is zero._{Y}**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`W`

or cross-identifiability gramian_{I}{2}`W`

are zero._{J}{2}**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.

### See Also

- Model Reduction Routines (Another Empirical Gramian Software; empirical W
_{C}and W_{O}only) - gram (Octave Control Package; linear W
_{C}and W_{O}only) - empirical gramian reference list

### Category

- MSC: 93A15, 93B11, 93C10
- ACM: Mathematical Software
- PACS: 02.30.Yy
- Controllability, Observability
- Science/Math/Software/MATLAB
- Programming Languages: Matlab
- math.OC, cs.SY, cs.MS
- Topic :: Scientific/Engineering :: Mathematics

### tl;dr

`emgr`

computes empirical system Gramian matrices that quantify *controllability*, *observability*, *minimality* or *identifiability* of nonlinear input-output systems.