# A Unified Software Framework for Empirical Gramians

Authors: Christian Himpe; Mario Ohlberger
Journal: Journal of Mathematics
 (2013)

---

## Abstract

A common approach in model reduction is balanced truncation, which is based on Gramian matrices classifying 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 provided 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 systems.


---

## Body

## 1. Introduction

In a control system setting, balanced truncation is a well-known 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 by
(1)x˙(t)=f(x,u,p),y(t)=g(x,u,p),
with the system function f and output function g of states x, input u, and parameters p; in the special case of an unparameterized linear system f=Ax(t)+Bu(t) and g=Cx(t). These empirical Gramians are computed by averaging simulations or experimental data with perturbations in inputs and initial states.Theemgr 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).


## 2. 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:(2)x˙(t)=Ax(t)+Bu(t),y(t)=Cx(t),
with the states x∈ℝn, control or input u∈ℝm, output y∈ℝo, system matrix A∈ℝn×n, input matrix B∈ℝn×m, and output matrix C∈ℝo×n.The necessary perturbations are given by six sets, of which{Eu,Ru,Qu} define the input perturbations, while sets {Ex,Rx,Qx} define the initial state perturbations:
(3)Eu={ei∈ℝj;∥ei∥=1;eiej≠i=0;i=1,…,m},Ex={fi∈ℝn;∥fi∥=1;fifj≠i=0;i=1,…,n},Ru={Si∈ℝj×j;Si*Si=1;i=1,…,s},Rx={Ti∈ℝn×n;Ti*Ti=1;i=1,…,t},Qu={ci∈ℝ;ci>0;i=1,…,q},Qx={di∈ℝ;di>0;i=1,…,r}.
These sets should correspond to the ranges in inputs and initial states the system is operating in.


### 2.1. 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:AWC+WCAT=-BBT. If the underlying system is asymptotically stable, the controllability Gramian can also be defined using the linear input-to-state map:
(4)WC=∫0∞eAτBBTeATτdτ.
Following [3], the empirical controllability Gramian is defined by the following.Definition 1.
For setsEu, Ru, and Qu, and input u(t), and input during the steady state x-(u-), the empirical controllability Gramianis given by
(5)WC=1|Qu||Ru|∑h=1|Qu|‍∑i=1|Ru|‍∑j=1m‍1ch2∫0∞Ψhij(t)dt,Ψhij(t)=(xhij(t)-x-)(xhij(t)-x-)*∈ℝn×n,
with xhij being the states for the input configuration uhij(t)=chSieju(t)+u-.Originally, in [2], u(t) was restricted to δ(t) but extended in [3] to arbitrary input under the name of empirical covariance matrix. x- can be the arithmetic average, the median, the steady state, or the principal components. Restricting Ru to {-1,1} simplifies the input perturbation to
(6)uhij(t)=-1iqheju(t)+u-.


### 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:ATWO+WOA=-CTC. Given an asymptotically stable underlying system, the observability Gramian can also be defined using the state-to-output map:
(7)WO=∫0∞eATτCTCeAτdτ.
The empirical observability Gramian is defined as described in [2, 3].Definition 2.
For setsEx, Rx, and Qx and output y during the steady state x-(y-), the empirical observability Gramianis given by
(8)WO=1|Qx||Rx|∑k=1|Qx|‍∑l=1|Rx|‍1dk2Tl∫0∞Ψkl(t)dtTl*,Ψabkl=(ykla(t)-y-)*(yklb(t)-y-)∈ℝ,
with ykla being the systems output for the initial state configuration x0kla=dkSlfa+x-.y
- can be the arithmetic average, the median, the steady-state output, or the principal components. Restricting Rx to {-1,1} simplifies (8) to
(9)WO=1|Qx||Rx|∑k=1|Qx|‍∑l=1|Rx|‍1dk2∫0∞Ψkl(t)dt,
and simplifies the initial state perturbation to
(10)x0kla=-1ldkfa+x-.


### 2.3. Cross Gramian

The cross Gramian [4] makes a combined 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: AWX+WXA=-BC. The cross Gramian can also be defined using the input-to-state and state-to-output maps, if the underlying system is asymptotically stable:
(11)WX=∫0∞eAτBCeAτdτ.
The empirical cross Gramian has been introduced in [5] for SISO systems and was extended to MIMO systems in [6].Definition 3.
For setsEu, Ex, Ru, Rx, Qu, and Qx and input u- during steady state x- with output y-, the empirical cross Gramian is given by
(12)WX=1|Qu||Ru|m|Qx||Rx|·∑h=1|Qu|‍∑i=1|Ru|‍∑j=1m‍∑k=1|Qx|‍∑l=1|Rx|‍1chdk∫0∞TlΨhijkl(t)Tl*dt,Ψabhijkl(t)=fb*Tk*Δxhij(t)ei*Sh*Δykla(t),Δxhij(t)=(xhij(t)-x-),Δykla(t)=(ykla(t)-y-),
with xhij and ykla being the states and output for the input uhij(t)=chSieju(t)+u- and initial state x0kla=dkTlfa+x-, respectively.x
- and y- can be the arithmetic average, the median, the steady state, or the principal components of the output. Again, restricting Ru and Rx to {-1,1} simplifies (12) to
(13)Ψabhijkl=(-1)i+lfb*Δxhij(t)ej*Δykla(t),
as well as simplifying input and initial state perturbation to
(14)uhij(t)=-1icheju(t)+u-,x0kla=-1ldkfa+x-.


### 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:
(15)x˙=f(x,u,p)=f(x,u)+∑k=1P‍f(x,pk).
The parameters p∈ℝP 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 WC,k. The sensitivity Gramian is now given by the diagonal matrix, with each diagonal element being the trace of a subcontrollability Gramian:
(16)WS=∑k=1P‍tr(WC,k)(δij)i=j=k.
The controllability and thus identifiability of each parameter are then given by the corresponding diagonal entry of the sensitivity Gramian WS. For partitionable linear systems, the sum of all subsystems controllability Gramians WC,k and the parameter-free subsystems Gramian WC,0 equals the usual controllability Gramian [7]:
(17)WC=WC,0+∑k=1P‍WC,k.
The sensitivity Gramian can be applied to nonpartitionable models with reduced accuracy.


### 2.5. 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 has(18)x˘˙=(x˙p˙)=(f(x,u,p)0),y=g(x,u,p).
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:
(19)W=(WOWQWQ*WP)⟹WI=WP-WQ*WO-1WQ.


### 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.
(20)x˘˙=(x˙p˙)=(f(x,u,p)v),y=(g(x,u,p)1),
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 WI¨ which can be extracted with a Schur complement from the joint Gramian:
(21)WJ=(WXWQWqWP)⟹WI¨=WP-WQ*WX-1WQ.


## 2.1. 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:AWC+WCAT=-BBT. If the underlying system is asymptotically stable, the controllability Gramian can also be defined using the linear input-to-state map:
(4)WC=∫0∞eAτBBTeATτdτ.
Following [3], the empirical controllability Gramian is defined by the following.Definition 1.
For setsEu, Ru, and Qu, and input u(t), and input during the steady state x-(u-), the empirical controllability Gramianis given by
(5)WC=1|Qu||Ru|∑h=1|Qu|‍∑i=1|Ru|‍∑j=1m‍1ch2∫0∞Ψhij(t)dt,Ψhij(t)=(xhij(t)-x-)(xhij(t)-x-)*∈ℝn×n,
with xhij being the states for the input configuration uhij(t)=chSieju(t)+u-.Originally, in [2], u(t) was restricted to δ(t) but extended in [3] to arbitrary input under the name of empirical covariance matrix. x- can be the arithmetic average, the median, the steady state, or the principal components. Restricting Ru to {-1,1} simplifies the input perturbation to
(6)uhij(t)=-1iqheju(t)+u-.


## 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:ATWO+WOA=-CTC. Given an asymptotically stable underlying system, the observability Gramian can also be defined using the state-to-output map:
(7)WO=∫0∞eATτCTCeAτdτ.
The empirical observability Gramian is defined as described in [2, 3].Definition 2.
For setsEx, Rx, and Qx and output y during the steady state x-(y-), the empirical observability Gramianis given by
(8)WO=1|Qx||Rx|∑k=1|Qx|‍∑l=1|Rx|‍1dk2Tl∫0∞Ψkl(t)dtTl*,Ψabkl=(ykla(t)-y-)*(yklb(t)-y-)∈ℝ,
with ykla being the systems output for the initial state configuration x0kla=dkSlfa+x-.y
- can be the arithmetic average, the median, the steady-state output, or the principal components. Restricting Rx to {-1,1} simplifies (8) to
(9)WO=1|Qx||Rx|∑k=1|Qx|‍∑l=1|Rx|‍1dk2∫0∞Ψkl(t)dt,
and simplifies the initial state perturbation to
(10)x0kla=-1ldkfa+x-.


## 2.3. Cross Gramian

The cross Gramian [4] makes a combined 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: AWX+WXA=-BC. The cross Gramian can also be defined using the input-to-state and state-to-output maps, if the underlying system is asymptotically stable:
(11)WX=∫0∞eAτBCeAτdτ.
The empirical cross Gramian has been introduced in [5] for SISO systems and was extended to MIMO systems in [6].Definition 3.
For setsEu, Ex, Ru, Rx, Qu, and Qx and input u- during steady state x- with output y-, the empirical cross Gramian is given by
(12)WX=1|Qu||Ru|m|Qx||Rx|·∑h=1|Qu|‍∑i=1|Ru|‍∑j=1m‍∑k=1|Qx|‍∑l=1|Rx|‍1chdk∫0∞TlΨhijkl(t)Tl*dt,Ψabhijkl(t)=fb*Tk*Δxhij(t)ei*Sh*Δykla(t),Δxhij(t)=(xhij(t)-x-),Δykla(t)=(ykla(t)-y-),
with xhij and ykla being the states and output for the input uhij(t)=chSieju(t)+u- and initial state x0kla=dkTlfa+x-, respectively.x
- and y- can be the arithmetic average, the median, the steady state, or the principal components of the output. Again, restricting Ru and Rx to {-1,1} simplifies (12) to
(13)Ψabhijkl=(-1)i+lfb*Δxhij(t)ej*Δykla(t),
as well as simplifying input and initial state perturbation to
(14)uhij(t)=-1icheju(t)+u-,x0kla=-1ldkfa+x-.


## 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:
(15)x˙=f(x,u,p)=f(x,u)+∑k=1P‍f(x,pk).
The parameters p∈ℝP 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 WC,k. The sensitivity Gramian is now given by the diagonal matrix, with each diagonal element being the trace of a subcontrollability Gramian:
(16)WS=∑k=1P‍tr(WC,k)(δij)i=j=k.
The controllability and thus identifiability of each parameter are then given by the corresponding diagonal entry of the sensitivity Gramian WS. For partitionable linear systems, the sum of all subsystems controllability Gramians WC,k and the parameter-free subsystems Gramian WC,0 equals the usual controllability Gramian [7]:
(17)WC=WC,0+∑k=1P‍WC,k.
The sensitivity Gramian can be applied to nonpartitionable models with reduced accuracy.


## 2.5. 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 has(18)x˘˙=(x˙p˙)=(f(x,u,p)0),y=g(x,u,p).
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:
(19)W=(WOWQWQ*WP)⟹WI=WP-WQ*WO-1WQ.


## 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.
(20)x˘˙=(x˙p˙)=(f(x,u,p)v),y=(g(x,u,p)1),
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 WI¨ which can be extracted with a Schur complement from the joint Gramian:
(21)WJ=(WXWQWqWP)⟹WI¨=WP-WQ*WX-1WQ.


## 3. Implementation

Theemgr software framework presented here provides a uniform interface to compute all six empirical Gramians and is given by
(22)W=emgr(f,g,q,t,w,p,vcfg,u,us,xs,um,xm,yd),
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.Table 1
Empirical gramian application matrix.
Gramian type
State reduction
Parameter reduction
Combined reduction
W
C
(✓)
×
×
W
O
(✓)
×
×
W
C
&
W
O
✓
×
×
W
X
✓
(✓)
×
W
S
(✓)
✓
×
W
I
(✓)
✓
(✓)
W
J
(✓)
✓
✓The following arguments are optional.p holds any parameters. The ten-component vector vcfg configures the available options, including averaging types, input and state scale subdivisions, and perturbation rotations. u provides the input to f and g, while us and xs set steady input and steady state. um and xm define the scales of the perturbation. Lastly, yd 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
(23)Ψabkl=(ykla(t)-y-)*(yklb(t)-y-)⟶{ψakl(t)=(ykla(t)-y-)Ψabkl=vec(ψakl)*vec(ψbkl).
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 at http://gramian.de or at the Oberwolfach References on Mathematical Software.


## 4. Numerical Experiments

To demonstrate the various empirical Gramians, computed by theemgr 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 B is also a random matrix, and the output matrix is given by C=BT. Furthermore a random, but element-wise, parametrized source term Ep of dimension n=100 parameters is added. Input is applied through a delta impulse:
(24)f(x,u,p)=x˙=Aarsinh(x)+Bu+Ep,g(x)=y=Cx.First, a state reduction, using the empirical controllability Gramian and the empirical observability Gramian, through balanced truncation is performed in Figure1, 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.Figure 1
Relative error in reduced system output by balanced truncation using the empirical controllability Gramian and the empirical observability GramianWC = emgr(f,g,[10,100,10],[0,0.01,1],“c”,p); WO = emgr(f,g,[10,100,10],[0,0.01,1],“o”,p).Next, a state reduction by direct truncation employing the empirical cross Gramian is demonstrated in Figure2, 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.Figure 2
Relative error in reduced system output by truncation using the empirical cross GramianWX = emgr(f,g,[10,100,10],[0,0.01,1],“x”,p).The empirical sensitivity Gramian can be applied if the underlying system can be partitioned such thatf(x,u,p)=f(x,u)+f(x,p). 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.Figure 3
Relative error in reduced system output with reduced source term; reduction by truncation using the empirical sensitivity GramianWS = emgr(f,g,[10,100,10],[0,0.01,1],“s”,p).Figure 4
Relative error in reduced system output with reduced source term; reduction by truncation using the empirical identifiability GramianWI = emgr(f,g,[10,100,10],[0,0.01,1],“i”,p).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 Figure5, 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 reduced-order model, of which the relative error is comparable to the other reduced system output models.Figure 5
Relative error in system output with combined state and parameter reduction; reduction by truncation of parameters and states using the empirical joint GramianWJ = emgr(f,g,[10,100,10],[0,0.01,1],“j”,p).


## 5. Future Work

Theemgr 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.



---

*Source: 365909-2013-09-19.xml*