A Unified Software Framework for Empirical Gramians

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.

Academic Editor:Beny Neta

Institute for Computational and Applied Mathematics, University of Münster, Einsteinstraße 62, 48149 Münster, Germany

Received 11 July 2013; Accepted 2 August 2013

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 [figure omitted; refer to PDF] with the system function [figure omitted; refer to PDF] and output function [figure omitted; refer to PDF] of states [figure omitted; refer to PDF] , input [figure omitted; refer to PDF] , and parameters [figure omitted; refer to PDF] ; in the special case of an unparameterized linear system [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . 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 [figure omitted; refer to PDF] (unit matrix) and [figure omitted; refer to PDF] (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: [figure omitted; refer to PDF] with the states [figure omitted; refer to PDF] , control or input [figure omitted; refer to PDF] , output [figure omitted; refer to PDF] , system matrix [figure omitted; refer to PDF] , input matrix [figure omitted; refer to PDF] , and output matrix [figure omitted; refer to PDF] .

The necessary perturbations are given by six sets, of which [figure omitted; refer to PDF] define the input perturbations, while sets [figure omitted; refer to PDF] define the initial state perturbations: [figure omitted; refer to PDF] 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: [figure omitted; refer to PDF] . If the underlying system is asymptotically stable, the controllability Gramian can also be defined using the linear input-to-state map: [figure omitted; refer to PDF] Following [3], the empirical controllability Gramian is defined by the following.

Definition 1.

For sets [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , and input [figure omitted; refer to PDF] , and input during the steady state [figure omitted; refer to PDF] , the empirical controllability Gramian is given by [figure omitted; refer to PDF] with [figure omitted; refer to PDF] being the states for the input configuration [figure omitted; refer to PDF] .

Originally, in [2], [figure omitted; refer to PDF] was restricted to [figure omitted; refer to PDF] but extended in [3] to arbitrary input under the name of empirical covariance matrix. [figure omitted; refer to PDF] can be the arithmetic average, the median, the steady state, or the principal components. Restricting [figure omitted; refer to PDF] to [figure omitted; refer to PDF] simplifies the input perturbation to [figure omitted; refer to PDF]

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: [figure omitted; refer to PDF] . Given an asymptotically stable underlying system, the observability Gramian can also be defined using the state-to-output map: [figure omitted; refer to PDF] The empirical observability Gramian is defined as described in [2, 3].

Definition 2.

For sets [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] and output [figure omitted; refer to PDF] during the steady state [figure omitted; refer to PDF] , the empirical observability Gramian is given by [figure omitted; refer to PDF] with [figure omitted; refer to PDF] being the systems output for the initial state configuration [figure omitted; refer to PDF] .

[figure omitted; refer to PDF] can be the arithmetic average, the median, the steady-state output, or the principal components. Restricting [figure omitted; refer to PDF] to [figure omitted; refer to PDF] simplifies (8) to [figure omitted; refer to PDF] and simplifies the initial state perturbation to [figure omitted; refer to PDF]

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: [figure omitted; refer to PDF] . The cross Gramian can also be defined using the input-to-state and state-to-output maps, if the underlying system is asymptotically stable: [figure omitted; refer to PDF] The empirical cross Gramian has been introduced in [5] for SISO systems and was extended to MIMO systems in [6].

Definition 3.

For sets [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] and input [figure omitted; refer to PDF] during steady state [figure omitted; refer to PDF] with output [figure omitted; refer to PDF] , the empirical cross Gramian is given by [figure omitted; refer to PDF] with [figure omitted; refer to PDF] and [figure omitted; refer to PDF] being the states and output for the input [figure omitted; refer to PDF] and initial state [figure omitted; refer to PDF] , respectively.

[figure omitted; refer to PDF] and [figure omitted; refer to PDF] can be the arithmetic average, the median, the steady state, or the principal components of the output. Again, restricting [figure omitted; refer to PDF] and [figure omitted; refer to PDF] to [figure omitted; refer to PDF] simplifies (12) to [figure omitted; refer to PDF] as well as simplifying input and initial state perturbation to [figure omitted; refer to PDF]

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: [figure omitted; refer to PDF] The parameters [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] . The sensitivity Gramian is now given by the diagonal matrix, with each diagonal element being the trace of a subcontrollability Gramian: [figure omitted; refer to PDF] The controllability and thus identifiability of each parameter are then given by the corresponding diagonal entry of the sensitivity Gramian [figure omitted; refer to PDF] . For partitionable linear systems, the sum of all subsystems controllability Gramians [figure omitted; refer to PDF] and the parameter-free subsystems Gramian [figure omitted; refer to PDF] equals the usual controllability Gramian [7]: [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] 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: [figure omitted; refer to PDF]

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. [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] which can be extracted with a Schur complement from the joint Gramian: [figure omitted; refer to PDF]

3. Implementation

The emgr software framework presented here provides a uniform interface to compute all six empirical Gramians and is given by [figure omitted; refer to PDF] 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
[figure omitted; refer to PDF]	([checkmark] )	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
[figure omitted; refer to PDF]	([checkmark] )	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
[figure omitted; refer to PDF]	[checkmark]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
[figure omitted; refer to PDF]	[checkmark]	([checkmark] )	[figure omitted; refer to PDF]
[figure omitted; refer to PDF]	([checkmark] )	[checkmark]	[figure omitted; refer to PDF]
[figure omitted; refer to PDF]	([checkmark] )	[checkmark]	([checkmark] )
[figure omitted; refer to PDF]	([checkmark] )	[checkmark]	[checkmark]

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 [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , 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 [figure omitted; refer to PDF] 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 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 [figure omitted; refer to PDF] is also a random matrix, and the output matrix is given by [figure omitted; refer to PDF] . Furthermore a random, but element-wise, parametrized source term [figure omitted; refer to PDF] of dimension [figure omitted; refer to PDF] parameters is added. Input is applied through a delta impulse: [figure omitted; refer to PDF]

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.

Figure 1: Relative error in reduced system output by balanced truncation using the empirical controllability Gramian and the empirical observability Gramian [figure omitted; refer to PDF] = emgr(f,g, [figure omitted; refer to PDF] 10,100,10 [figure omitted; refer to PDF] , [figure omitted; refer to PDF] 0,0.01,1 [figure omitted; refer to PDF] ,"c",p); W_O = emgr(f,g, [figure omitted; refer to PDF] 10,100,10 [figure omitted; refer to PDF] , [figure omitted; refer to PDF] 0,0.01,1 [figure omitted; refer to PDF] ,"o",p).

[figure omitted; refer to PDF]

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.

Figure 2: Relative error in reduced system output by truncation using the empirical cross Gramian [figure omitted; refer to PDF] = emgr(f,g, [figure omitted; refer to PDF] 10,100,10 [figure omitted; refer to PDF] , [figure omitted; refer to PDF] 0,0.01,1 [figure omitted; refer to PDF] ,"x",p).

[figure omitted; refer to PDF]

The empirical sensitivity Gramian can be applied if the underlying system can be partitioned such that [figure omitted; refer to PDF] . 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 Gramian [figure omitted; refer to PDF] = emgr(f,g, [figure omitted; refer to PDF] 10,100,10 [figure omitted; refer to PDF] , [figure omitted; refer to PDF] 0,0.01,1 [figure omitted; refer to PDF] ,"s",p).

[figure omitted; refer to PDF]

Figure 4: Relative error in reduced system output with reduced source term; reduction by truncation using the empirical identifiability Gramian [figure omitted; refer to PDF] = emgr(f,g, [figure omitted; refer to PDF] 10,100,10 [figure omitted; refer to PDF] , [figure omitted; refer to PDF] 0,0.01,1 [figure omitted; refer to PDF] ,"i",p).

[figure omitted; refer to PDF]

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 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 Gramian [figure omitted; refer to PDF] = emgr(f,g, [figure omitted; refer to PDF] 10,100,10 [figure omitted; refer to PDF] , [figure omitted; refer to PDF] 0,0.01,1 [figure omitted; refer to PDF] ,"j",p).

[figure omitted; refer to PDF]

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

Acknowledgments

The authors acknowledge the support by the Deutsche Forschungsgemeinschaft and the Open Access Publication Fund of the University of Münster.

References

[1] B. C. Moore, "Principal component analysis in linear systems: controllability, observability, and model reduction," IEEE Transactions on Automatic Control , vol. 26, no. 1, pp. 17-32, 1981.

[2] S. Lall, J. E. Marsden, S. Glavaski, "Empirical model reduction of controlled nonlinear systems," in Proceedings of the International Federation of Automatic Control World Congress (IFAC '99), pp. 473-478, 1999.

[3] J. Hahn, T. F. Edgar, "Balancing approach to minimal realization and model reduction of stable nonlinear systems," Industrial and Engineering Chemistry Research , vol. 41, no. 9, pp. 2204-2212, 2002.

[4] A. C. Antoulas Approximation of Large-Scale Dynamical Systems , of Advances in Design and Control, pp. xxvi+479, SIAM, Philadelphia, Pa, USA, 2005.

[5] S. Streif, R. Findeisen, E. Bullinger, "Relating cross gramians and sensitivity analysis in systems biology," in Proceedings of the 17th International Symposium on Mathematical Theory of Networks and Systems (MTNS '06), pp. 437-442, 2006.

[6] C. Himpe, M. Ohlberger, "Cross-Gramian based combined state and parameter reduction," http://arxiv.org/abs/1302.0634

[7] C. Sun, J. Hahn, "Parameter reduction for stable dynamical systems based on Hankel singular values and sensitivity analysis," Chemical Engineering Science , vol. 61, no. 16, pp. 5393-5403, 2006.

[8] D. Geffen, R. Findeisen, M. Schliemann, F. Allgöwer, M. Guay, "Observability based parameter identifiability for biochemical reaction networks," in Proceedings of the American Control Conference (ACC '08), pp. 2130-2135, June 2008.