A Unified Software Framework for Empirical Gramians

A common approach in model reduction is balanced truncation, which is based on gramian matrices classifiying certain attributes of states or parameters of a given dynamic system. Initially restricted to linear systems, the empirical gramians not only extended this concept to nonlinear systems, but also provide a uniform computational method. This work introduces a unified software framework supplying routines for six types of empirical gramians. The gramian types will be discussed and applied in a model reduction framework for multiple-input-multiple-output (MIMO) systems.


Introduction
In a control system setting, balanced truncation is a wellknown technique for model reduction. Introduced by [1], Gramian matrices were employed to determine controllability and observability of linear systems. From these Gramians, a balancing transformation can be computed, enabling the truncation, for example, of states that are neither controllable nor observable.
With [2], empirical (controllability and observability) Gramians were introduced, which correspond to the analytical Gramians for linear systems, while extending the concept of system Gramians to nonlinear systems which are generally given bẏ( with the system function and output function of states , input , and parameters ; in the special case of an unparameterized linear system = ( ) + ( ) and = ( ). These empirical Gramians are computed by averaging simulations or experimental data with perturbations in inputs and initial states.
The emgr framework presented here encompasses six empirical Gramians, namely, the controllability, observability, cross, sensitivity, identifiability, and joint Gramians. To adapt the computation of empirical Gramians to the operating setting of the system, the initial state and the input are the main variables which are perturbed by rotations and scaling. The sets of rotations provided are {1} (unit matrix) and {−1, 1} (negative unit matrix and unit matrix). Though these are very basic sets and thus might not reflect all dynamics, especially with interrelated states and parameters, they allow a very efficient Gramian assembly. Scales may be freely chosen. The subdivision of the scales may be linear, logarithmic, or geometric. Finally, there are several options to average against the arithmetic average [2], the median, a steady state [3], and additionally, the principal components of the simulations or data via a proper orthogonal decomposition (POD).

Empirical Gramians
Concerned with the reduction of states, the controllability, observability, and cross Gramians are presented next, followed by the sensitivity, identifiability, and joint Gramians, which are used for parameter and combined reduction. For the purpose of defining the Gramians, a linear time-invariant control system is assumed: with the states ∈ R , control or input ∈ R , output ∈ R , system matrix ∈ R × , input matrix ∈ R × , and output matrix ∈ R × .
The necessary perturbations are given by six sets, of which { , , } define the input perturbations, while sets { , , } define the initial state perturbations: ( These sets should correspond to the ranges in inputs and initial states the system is operating in.

Controllability Gramian.
Controllability is a quantification of how well a state can be driven by input. Analytically, the controllability Gramian is given by the smallest semipositive definite solution of the Lyapunov equation: . If the underlying system is asymptotically stable, the controllability Gramian can also be defined using the linear input-to-state map: (4) Following [3], the empirical controllability Gramian is defined by the following.
Definition 1. For sets , , and , and input ( ), and input during the steady state ( ), the empirical controllability Gramian is given by with ℎ being the states for the input configuration ℎ ( ) = ℎ ( ) + .
Originally, in [2], ( ) was restricted to ( ) but extended in [3] to arbitrary input under the name of empirical covariance matrix. can be the arithmetic average, the median, the steady state, or the principal components. Restricting to {−1, 1} simplifies the input perturbation to 2.2. Observability Gramian. Observability quantifies how well a change in a state is reflected by the output. The analytical observability Gramian is given by the smallest semipositive definite solution of the Lyapunov equation: Given an asymptotically stable underlying system, the observability Gramian can also be defined using the state-to-output map: The empirical observability Gramian is defined as described in [2,3].
Definition 2. For sets , , and and output during the steady state ( ), the empirical observability Gramian is given by with being the systems output for the initial state configuration 0 = + .
can be the arithmetic average, the median, the steadystate output, or the principal components. Restricting to and simplifies the initial state perturbation to

Cross Gramian. The cross Gramian [4] makes a com-
bined statement about the controllability and observability, given that the system has the same number of inputs and outputs. If the system is also symmetric, meaning that the system transfer function is symmetric, then the absolute value of these Gramians' eigenvalues equals the Hankel singular values. It is originally computed as the smallest solution of the Sylvester equation: + = − . The cross Gramian can also be defined using the input-to-state and state-tooutput maps, if the underlying system is asymptotically stable: The empirical cross Gramian has been introduced in [5] for SISO systems and was extended to MIMO systems in [6].
Definition 3. For sets , , , , , and and input during steady state with output , the empirical cross Gramian is given by with ℎ and being the states and output for the input ℎ ( ) = ℎ ( ) + and initial state 0 = + , respectively. and can be the arithmetic average, the median, the steady state, or the principal components of the output. Again, restricting and to {−1, 1} simplifies (12) to as well as simplifying input and initial state perturbation to 2.4. Sensitivity Gramian. The sensitivity Gramian allows controllability-based parameter reduction and identification.
It is based on [7] and aimed for models that can be partitioned as follows:̇= The parameters ∈ R are handled here as additional inputs. All summands of the partitioned system are treated as independent subsystems, and thus a controllability Gramian for each subsystem can be computed. Each parameter controllability is encoded in the sum of singular values of the associated subcontrollability Gramian , . The sensitivity Gramian is now given by the diagonal matrix, with each diagonal element being the trace of a subcontrollability Gramian: The controllability and thus identifiability of each parameter are then given by the corresponding diagonal entry of the sensitivity Gramian . For partitionable linear systems, the sum of all subsystems controllability Gramians , and the parameter-free subsystems Gramian ,0 equals the usual controllability Gramian [7]: The sensitivity Gramian can be applied to nonpartitionable models with reduced accuracy.

Identifiability
Gramian. The identifiability Gramian enables observability-based parameter identification and consequently parameter reduction. As described in [8], the dynamic system states are augmented with as many states as parameters that are constant over time and have the initial value of the (prior) parameter value. One haṡ̆= The observability Gramian of this augmented system holds the observability information of states and parameters. To extract the parameter specific observability, the Schur complement can be applied to the augmented observability Gramian: 2.6. Joint Gramian. Based on the identifiability Gramian procedure, the cross Gramian can be employed for a concurrent state and parameter reduction (see [6]). As for the identifiability gramian, the system is augmented by constant parameter states. Additionally, as many inputs and outputs as parameters are augmented as well. The augmented inputs act via identity on the augmented states. Likewise, the augmented states are mapped by identity to the augmented outputs to preserve symmetry.̆= to preserve symmetry. The cross Gramian of this special augmented system, similar to the identifiability Gramian, holds the cross Gramian of the original system as well as a cross-identifiability Gramian̈which can be extracted with a Schur complement from the joint Gramian:

Implementation
The emgr software framework presented here provides a uniform interface to compute all six empirical Gramians and is given by with f and g being handles to the system and the output function, both requiring the signature f(x,u,p) and g(x,u,p). q is a vector defining the systems number of inputs, states, and outputs. t is a three-component vector containing start time, time step, and stop time. w is a character setting the Gramian type; for an overview on the applicability of Gramian types, see Table 1.
The following arguments are optional. p holds any parameters. The ten-component vector v cfg configures the available options, including averaging types, input and state scale subdivisions, and perturbation rotations. u provides the input to and , while u s and x s set steady input and steady state. u m and x m define the scales of the perturbation. Lastly, y d allows passing experimental data to be used instead of generated snapshots.
The parameters reducing empirical Gramians (sensitivity, identifiability, and joint) are an encapsulation of the state reducing empirical Gramians (controllability, observability, and cross). Computation of the latter is extensively vectorized, exploiting the Gramian matrix assembly format. For example, the empirical observability Gramian assembly from (8) can be computationally simplified to Computation of empirical Gramians using emgr is very portable, since only basic vector and matrix operations are required. Necessary integrations, meaning simulations for given inputs or initial states, are accomplished either by the first-order Euler's method, second-order Adams-Bashforth method, or second-order Leapfrog method. The empirical Gramian framework emgr as well as the following experiments is released under an open source license, is compatible with OCTAVE and MATLAB, and can be found

Numerical Experiments
To demonstrate the various empirical Gramians, computed by the emgr framework, a symmetric nonlinear MIMO system with one hundred states, ten inputs, and ten outputs is employed. The system matrix is generated randomly with ensured stability and symmetry; the input matrix is also a random matrix, and the output matrix is given by = . Furthermore a random, but element-wise, parametrized source term of dimension = 100 parameters is added. Input is applied through a delta impulse: (24) First, a state reduction, using the empirical controllability Gramian and the empirical observability Gramian, through balanced truncation is performed in Figure 1, reducing the number of states to the number of outputs. Balanced truncation as a classic approach in model order reduction will be used as a baseline to which the following methods will be compared.
Next, a state reduction by direct truncation employing the empirical cross Gramian is demonstrated in Figure 2, again reducing the number of states to ten. The state reduction via direct truncation of the cross Gramian has about the same error but requires only one Gramian and no balancing.
The empirical sensitivity Gramian can be applied if the underlying system can be partitioned such that ( , , ) = ( , ) + ( , ). To be able to use it in this setting, the parametrized source term is reduced to the number of outputs in Figure 3. The sensitivity Gramian is the fastest parameter reduction method but has a high relative error in outputs. Since the parameters of the source term are reduced, the cumulative effects in the original system are the origin of the increasing error over time. Next, the parametrized source term is reduced by the empirical identifiability Gramian in Figure 4.
Taking five times as long for the parameter reduction as the sensitivity Gramian, the identifiability Gramian is about two orders of magnitude more accurate.
Finally, in Figure 5, the same system undergoes a combined state and parameter reduction using the empirical joint Gramian. The joint Gramian is the only Gramian allowing direct, balancing-free combined reduction of state and parameter space with a comparable relative error and yet takes about the same time for the reduction as the identifiability Gramian. This combined reduction generates a reducedorder model, of which the relative error is comparable to the other reduced system output models.

Future Work
The emgr framework already allows a wide range of computations of empirical Gramians for state or parameter reduction. Apart from model order reduction, the empirical Gramians can be employed for system identification tasks, like parameter identification or sensitivity analysis as well as decentralized control, nonlinearity measurement, and uncertainty quantification. Further work will enhance the flexibility, while keeping the interface as simple as possible. Following [8], allowing factorial designs will greatly enlarge the field of application. Finally, extending the use of the cross Gramian (and thus the joint Gramian) to nonsymmetric systems [4] will enable a combined state and parameter reduction for general linear and nonlinear models without balancing.