Order reduction of hierarchical interconnected dynamical systems

The simulation of large scale nonlinear dynamical interconnected systems, as they arise in all modern engineering disciplines, is a usual task. Due to the high complexity of the considered systems, the principle of thinking in hierarchical structures is essential and common among engineers. Therefore, this contribution proposes an approach for the numerical simulation of large systems, which keeps the hierarchical system structure alive during the entire simulation process while simultaneously exploiting it for order reduction purposes. This is accomplished by embedding the trajectory piecewise linear order reduction scheme in a modified variant of the component connection modeling for building interconnected system structures. The application of this concept is demonstrated by means of a widely used benchmark example and a modified variant of it.


Introduction
The design and operation of todays electrical systems, ranging from microelectronic circuits to power generation and distribution systems, is enabled by powerful simulation tools at hand.Models for such large scale systems are commonly built up hierarchically from smaller and less complex interconnected sub-systems.While this hierarchical point of view is the state-of-the-art, which is utilized within the modeling process itself, the actual simulation is performed by numerically integrating the resulting large scale system of nonlinear ordinary differential equations.
The same is true for model order reduction methods for linear systems, (cf.Antoulas et al., 2001;Antoulas, 2005;Baur et al., 2014), as well as for nonlinear systems, (cf.Baur et al., 2014;Rewienski and White, 2003a).All these methods start at a given description in form of an ordinary differential equation of the entire system to be reduced, dismissing possible hierarchical structures which could be exploited ben-eficially, e.g. by parallelizing the evaluation of the ordinary differential equation during a numerical integration for transient simulation.
In contrast to this common practice, the simulation approach proposed in this contribution makes use of the so called component connection modeling (CCM) by DeCarlo and Saeks (1981), which preserves the system's hierarchical structure as well as its internal interconnections during the entire modeling and simulation process.Details on the main conception of CCM and the actually applied modified version called mCCM (Popp et al., 2017a) with several structural changes to the original concept is introduced and discussed with respect to its properties in Sect. 4.
Besides the mCCM itself, this contribution demonstrates, how the trajectory piecewise-linear approach (TPWL) (Rewienski and White, 2003b) as an example for a broad variety of order reduction methods for nonlinear dynamical systems can be embedded into such a framework.The main benefit arising from this conceptual combination of mCCM and TPWL can be seen in the additional degree of freedom which is gained by exploiting the intrinsic hierarchy of large interconnected systems, i.e. the choice of the hierarchical level which is to be reduced in its order.This important degree of freedom refers to the necessity to obtain reference trajectories of the non-reduced original system as a starting point of the reduction process.Increasing the size of the considered systems, this necessity limits the practical applicability and the proposed hierarchical approach helps to overcome or at least to ease this problem.Therefore, Sect. 5 outlines the implemented concepts of CoSimMA, a fully functional simulation framework based on the mCCM formulation of interconnected dynamical systems.As a part of these concepts, the structural embedding of order reduction methods into the framework CoSimMA is considered.Two applicable procedures for reducing a model's order are briefly described and discussed with respect to their properties, advantages and drawbacks.The well established balanced truncation method is treated in Sect. 2 and for nonlinear systems the trajectory piecewise linear (TPWL) reduction scheme is chosen for further analysis in Sect.3, without restricting the general embedding concept in this work to these order reduction methods.
In Sect.6 the proposed combined approach of TPWL order reduction of a hierarchically built up model by using mCCM within the simulation framework CoSimMA is demonstrated.As an illustrative example, a standard benchmark system originally introduced in Chen and White (2000) is processed in the context of the mentioned methods.Furthermore, a slightly altered version of this example is utilized to show a significant sensitivity of a systems reducibility to the choice of output variables, which should be preserved during the order reduction process.

Order reduction of linear systems via balanced truncation
The order reduction of linear dynamical systems of the form with the order n, i.e. x = [x 1 , . .., x n ] , aims at finding an approximate system representation ẋr = A r x r + B r u (2) y r = C r x r of order q with x r = x r1 , . .., x rq and q n while approximately keeping the outputs unchanged, i. e. y ≈ y r .This order reduction is formally achieved by applying a state transformation wherein the transformation matrix T has the dimension q ×n and is assumed to be orthonormal, i. e. T T = 1 holds.Over the last decades several specialized methods for computing such an order reducing transformation in Eq. ( 3) have been developed and the properties of each method are well understood (cf.Antoulas et al., 2001;Antoulas, 2005;Baur et al., 2014 among many others).
As an example of such a method of computing a transformation matrix T, the balanced truncation (BT), often also referred to as truncated balanced realization (TBR), tracing back to Moore (1981) is discussed in the following.The BT is based on the controllability and the observability of the states x of the original system in Eq. (1), measured by the gramian matrices P and Q.According to Moore (1981), first a so called balancing transformation x b = T b x has to be found, which makes each state within x b itself equally well controllable and observable, which is expressed by wherein P b and Q b denote the gramian matrices of the balanced system, which are equal diagonal matrices with the Hankel singular values (HSV) σ 1 , . .., σ n in descending order.Algorithms for computing such a balancing transformation are given e.g. in Laub et al. (1987).States in x b with corresponding small HSV are neither well controllable nor well observable and can therefore be truncated in the second step of the BT procedure, resulting in a reduced order system in Eq. ( 2) of order q.
A key benefit of the BT reduction scheme described above is the existence of a maximum error bound (Glover, 1984) referred to the transfer functions H(s) and H r (s) of the original and the reduced system, respectively, which can be computed a priori, before setting up the actual reduced order system.With this at hand, the decision of how to choose the reduced order q can be made upon a quantitative criterion instead of a trial and error procedure, as it is the case with other reduction schemes (Antoulas, 2005).
3 Order reduction of nonlinear systems via TPWL For the reduction of nonlinear dynamical systems of the form the two basic steps of reducing linear systems, cf.Sect.2, can in principle be applied.Unfortunately, the computation of an appropriate state variable transformation as well as the reduction procedure itself are more challenging in the nonlinear case, strongly depending on the nonlinearity of the system in Eq. ( 6).
Restricting the structure of the considered original nonlinear dynamical system of order n to ẋ = g(x) + B(x)u (7) y = Cx and assuming the transformation x r = Tx as introduced in Sect. 2 to be appropriate and known, the resulting reduced nonlinear model of order q can be written as ẋr = Tg(T x r ) + TB(T x r )u (8) y r = C r x r with C r = CT .
While in the linear case, the system matrix of the reduced system of true order q is given by TAT , this does not work in the nonlinear case.Even if the terms Tg(T x r ) and TB(T x r ) in Eq. ( 8) appear to be of reduced order q, the computational complexity is even increased, because the nonlinear function g(•) and the state dependent input matrix B(•) and their arguments still have the original order n.The state transformation in general cannot be applied once during the reduction process.Instead, it has to be evaluated for each evaluation of the system's equations and thus no computational effort is saved in spite of the reduced model order (Rewienski and White, 2003b).
Over the last decades, several approaches have been developed in order to find usable state transformations T on the one hand and practical ways of avoiding the problem of an only seemingly reduced model order.For an overview of these developments (cf.Baur et al., 2014;Dukic and Saric, 2012) as comprehensive survey articles.
One of these approaches is the TPWL order reduction scheme originally developed by Rewienski and White (2003a), which underwent several extensions and variations (Dong and Roychowdhury, 2008;Vasilyev et al., 2003;Martinez, 2009).In this paper, only the original version Rewienski and White (2003a) with application of BT as the reduction method for the resulting linear systems similar to Vasilyev et al. (2003) is discussed.
In order to be able to evaluate the reduction transformations in Tg(T x r ) and TB(T x r ) for lowering the computational complexity as well as the model order, first, multiple linearizations of Eq. ( 7) around N different points x i in the state space are generated and afterwards these individual linearized models are put together in form of a convex combination which is called trajectory piecewise linear approximation of 7. The N state dependent functions w i (x) are weights, which select the most relevant linear models among the given set and holds.A i represents the Jacobian of g(x) evaluated at the linearization points x i and B i are the corresponding input matrices.
Since Eq. ( 9) consists of a collection of linear models, an order reducing state transformation x r = Tx can be applied and evaluated, resulting in as the reduced order TPWL model.Here, it is assumed, that the freedom of choosing the transformed weight functions can be used to render them evaluable with reasonable effort, to prevent the aforementioned problem associated with the order reduction of nonlinear systems.In Rewienski and White (2003a) such a choice for the weight function is given as which also can be applied as w ri (x r ) by substituting x r for x and Tx i for x i in Eq. ( 13).With an empirical choice of the constant parameter β the number of simultaneously active linear systems can be influenced.Also on an empirical basis, a procedure to select the linearization points x i based on euclidean distances in the state space can be found in Rewienski and White (2003a).In this approach, the full order system in Eq. 7 is simulated and each time the inequality min i=1,...,N is fulfilled, a new linearization point x i+1 is added.With the constant parameter δ, the final number N of linearization points can be controlled empirically, also strongly depending on the choice of the test trajectory.
Even if there are many examples in the literature, in which the TPWL (Rewienski and White, 2003a;Dong and Roychowdhury, 2008;Vasilyev et al., 2003;Martinez, 2009) and other order reduction schemes like Empirical Gramian based balanced trucation (Condon and Rossen, 2004) or Krylov subspace based reduction of quadratic taylor expanded nonlinear systems (Chen and White, 2000) among many others are applied successfully, i. e. the example system's order as well as the simulation time decreases notably, all methods have many empirical parameters and options.

Modified component connection modeling (mCCM)
Originally introduced by DeCarlo and Saeks (1981), CCM provides a mathematical framework to describe interconnected dynamical systems in a convenient way.It focuses on a strictly separate treatment of the individual sub-systems, called components, and the connections between them, which are used to build up an interconnected system.This separation of components and connections is kept during the entire modeling and simulation process.A comparison of CCM and other formalisms for the description of interconnected dynamical systems is given in Popp et al. (2017a).
Within the proposed mCCM, which is a structurally simplified variant of the original CCM aiming towards hierarchical system formulations, the interconnection process of individual components S i is performed in three stages and, optionally, hierarchical connections can be considered in a fourth stage.
In the first stage, the components S i , which are the smallest entities used for building up larger systems, have to be defined.This can be done by means of systems of explicit ordinary differential equations for dynamical components and as systems of algebraic equations for non-dnyamical systems, wherein u i is the vector of the input quantities and x i denotes the state vector of dynamical components in Eq. ( 15) or the output vector of nondynamical components in Eq. ( 16), respectively.Deviating from the original CCM, in mCCM more general descriptions of dynamical components like a mass matrix formulation M i (x i ) ẋi = f mi (x i , u i ) or even a fully implicit system of ordinary differential equations 0 = F i (x i , ẋi , u i ) are permitted, which both allow for the important case of differential algebraic component equations.For the sake of simplicity, only dynamical components in the form of Eq. ( 15) are considered in the following without loss of generality.As a further deviation from the original framework, any kind of observer equation is dismissed in the component description.This simplifies the formulation of connections between components as will be seen later in this section.
In the second formal step, the individual components S i , which are to be connected in order to form a larger system, are collated in the composite component S as depicted in Fig. 1a.Thereby, the composite component is described by collating the individual component's states, inputs and vector fields.Up to this step, no connections are considered, yet.
These are introduced in the third step by means of a linear algebraic connection equation in which the weighted connection between the i-th input variable in u and the j -th output variable in x is established by a nonzero matrix entry c ij within the connection matrix C.
The connected component C arises by combining the compound component S and the connection equation as illustrated in Fig. 1b.Besides the main properties from a system and network theoretic point of view, which are mentioned in DeCarlo and Saeks (1981) and Mathis (1987, pp.175 et seq.), all component input and output quantities are preserved throughout the entire modeling and simulation procedure.This is in contrast to the loss of intermediate information, which occurs by simply substituting component inputs by the corresponding linear combinations of outputs.Moreover, by dismissing observer equations, as mentioned above, a single connection between two components is expressed only by a single entry in the connection matrix C, which makes the formulation of hierarchically interconnected systems a straight forward procedure.
As depicted in Fig. 1c, not only individual components S i but also several connected components C i can be connected again, forming arbitrarily deep hierarchical model structures without any theoretical limitations.The hierarchical connected component C is formed by the same steps as described above.Only the formation of the hierarchical connection matrix C has to be handled additionally in a fourth step.This is accomplished by placing the connection matrices C i of the individual connected components C i as diagonal blocks in hierarchical connection matrix C and connections between individual connected components C i are handled as entries besides these diagonal blocks.At this point, it becomes clear how the structural simplifications of mCCM to the original concept in DeCarlo and Saeks (1981) make it very convenient, to model hierarchical system structures while keeping the key benefits of the CCM.

CoSimMA -an mCCM based simulation framework
As originally introduced in Popp et al. (2016a), CoSimMA, a modular simulation framework, uses mCCM as theoretical core for performing studies on the efficient simulation of large interconnected dynamical systems.The basic concept behind CoSimMA and its modular structure is graphically outlined in Fig. 2.After setting up the connected component, several types of model manipulation and simulation such as linearization, model aggregation, linear and nonlinear order reduction as well as static and transient simulations can be performed in a modular manner.The key of the implementation concept is, that the result of each manipulation is a connected component or a component, respectively.This is enabled by strictly utilizing the object oriented programming paradigm within the problem solving environment MATLAB (The MathWorks, Inc., 2015).Using MAT-LAB's interfaces for other programming languages such as C or FORTRAN allows the convenient usage of sophisticated software packages for simulation purposes such as SUN-DIALS (Hindmarsh et al., 2005) as an alternative to MAT-LAB's own numerical solvers for ordinary differential equations (Shampine et al., 2003;Shampine, 2002).On the one hand, due to the unified module interfaces within CoSimMA, it becomes completely transparent from a users point of view, which specific algorithm and software package is used for simulation purposes.On the other hand, the development of new simulator modules, such as the model order reduction methods treated in this paper, are easily implemented, since already existing modules are interfaced in an abstract and therefore very clear manner.All simulation results presented in the following section have been obtained with CoSimMA.

Nonlinear modeling and reduction example
In this section, first the application of the hierarchical modeling approach utilizing the mCCM framework is demonstrated by means of an example circuit.Afterwards, several reduced order models are generated by applying the TPWL order reduction scheme combined with the BT for reducing the intermediate linear systems.
The example circuit treated in this section originates from Chen and White (2000) and is used by several other authors, e.g.Rewienski and White (2003a), as a benchmark system for the comparison of different order reduction methods.It is composed of a current driven input stage followed by N nonlinear subcircuits in a chain configuration as depicted in Fig. 3a.Although this circuit is intended as an academic benchmark example, it can be used as a charge pump for the transformation of a low input voltage u C0 into a higher output voltage u CN .
In the original version of the considered example, the capacitor voltage u C0 of the input stage, cf.Fig. 3a, is chosen as output variable which is preserved throughout the reduc-  tion process.This specific choice apparently leads to a good reducibility of the complete system because chain sections located far away from this circuit node of interest are less important for the actual behavior of the chosen output quantity.
In order to demonstrate the severe influence of the choice of variables which are to be preserved throughout the reduction process, a second scenario is treated in the following.Therein, additionally to the fist capacitor voltage u C0 the voltage of the last capacitor in the chain network is chosen to be an output variable, i.e. y = [u C0 , u C50 ] for N = 50 as size of the original system.

Hierarchical modeling of the nonlinear example circuit
Following the presented steps of the mCCM modeling procedure in Sect.4, fist the components itself are defined.The input stage of the given circuit in Fig. 3a leads to a dynamical component C 0 depicted in Fig. 3b, which is described by the ordinary differential equation with the circuit element parameters R = 1 and C = 1F .All diodes of the entire circuit are described by the algebraic equation wherein k D = 40V −1 and I S = 1A are chosen.The vector of input variables of C 0 is given by u C0 = [i 0 , i e ] and the corresponding state vector results in x C0 = u C0 .The n-th chain section C n , shown in Fig. 3c, consists of a dynamical component C 1,n called chain element and a nondynamical observer component C 2,n called current observer.The component equation of the C 1,n reads with u C1,n = [u n , i n ] and x C1,n = u Cn as input and state vectors, respectively.The current input i n−1 of the chain section C n−1 is observed by means of the observer component C 2,n , which is described by and the corresponding input and state vectors are defined as u C2,n = [u n , u Cn ] and x C1,n = i n−1 .Both build up a connected component C n by applying a suitable connection equation to set up the interconnection structure within C n as shown in Fig. 3c.In order to compile the entire example system, a hierarchical connected component has to be set up by interconnecting the input element C 0 with N = 50 chain sections as depicted in Fig. 3d and adding the current source i q for the excitation of the system.
6.2 Reduced TPWL models 6.2.1 Reference simulation of the original system As described in Sect. 3 it is necessary to simulate the original system of full order n in order to obtain a training trajectory for setting up a TPWL reduced model.This reference simulation furthermore is used as a base for the normalized representation of computation times, i.e. the time needed to perform the reference simulation (210 ms) is referred to as 100 %.Even though the figures below show different time axis, all simulations have been performed for a time interval of 5 s with the same solver (ode45 from Shampine et al., 2003) for each resulting systems of ordinary differential equations.All relative errors stated below are also referred to the reference simulation in the entire considered time interval.

Reduced system -case 1
In this first reduction example the capacitor voltage u C0 of the input stage is chosen as the output variable, which is to be preserved after the reduction has been applied.This choice equals the original version of this example (Chen and White, 2000).
For the test trajectory resulting form zero initial conditions for all state variables and a step excitation of 1 A, 27 linearization points have been found empirically to be a good compromise between approximation quality and computational effort.The unreduced TPWL representation shows a maximal relative error smaller than 0.1 % and it increases the computation time by 35.2 % referred to the reference simulation.The behavior of the first five weight functions is shown in Fig. 4 and it can be seen, that most of the time only a single linearized model is active, as desired.
Reducing the order to q = 16 results in a maximal relative error smaller than 1% accompanied by a moderately decreased computation time by −9.9 %.A more significant reduction if the computation time by −61% can be reached with a reduced order of q = 1 while keeping the relative error below 3 %.Clearly, the usable range around the training trajectory of this reduced model is very limited.A comparison of the two reduced models and the original system is given in Fig. 5, revealing slight transient and static deviations within the stated a posteriori error bounds.

Reduced system -case 2
As a second case, a slightly modified variant of the original system in Chen and White (2000) is investigated.Instead of only preserving the capacitor voltage u C0 , in this case additionally u C50 is considered as an output variable.For the sake of comparison, the reduced order q had been lowered step by step in order to keep the relative error just below 1%.This could be reached with q = 34, which more then doubles the number or states in the reduced model compared to case 1.The simulation time increases by 46 %, which makes the order reduction superfluous for the practical use in this case.
The further reduction with q = 8 results in the first reduced model with actual savings in computation time (3.8 %), whereby the relative error increases to nearly 10 %.Comparing the small benefit in computation time with the increased relative error and the added slightly oscillatory behavior visible in the simulation results shown in Fig. 6 also reveals limited practical use of such a reduced model.
In order to clarify the severe differences in the reduction results between case 1 and case 2, the Hankel singular values of the linearized and balanced models of both cases are compared in Fig. 7.It turns out, that the Hankel singular values in case 1 decay very fast while in case 2 this decay is less steep.Recalling the relation between the Hankel singular values and the controllability and the observability of state variables explains, why the choice of u C50 as additional output variable lowers the reducibility of the entire system.

Conclusions
With the introduction of mCCM, a modified and structurally extended version of the component connection modeling by DeCarlo and Saeks, a useful framework for the convenient and efficient formulation of hierarchical interconnected dynamical systems had been presented and its advantageous properties were discussed in detail.Due to the rising computational complexity with an increasing number of components to be simulated, it is desirable to have advanced simulation techniques such as model order reduction methods in a usable form at hand.Therefore, it had been described, how such order reduction schemes can be embedded into the software framework CoSimMA, a fully functional simulator based on the mCCM formulation of interconnected dynamical systems.Two model order reduction schemes, namely the well established balanced truncation method for linear systems as well as the TPWL method for nonlinear systems were described and discussed with respect to their basic properties, advantages and disadvantages.Finally, a modeling, reduction and simulation study by means of a standard benchmark example had been presented and critically discussed with respect to the quality of the reduced models and the savings in computation time.
As a main result it had been shown, that the reducibility of a given system strongly depends on the choice of the output variables which are preserved throughout the reduction process.In the case of a well reducible system, considerable savings in computation time could be achieved while keeping the approximation error in moderate bounds.In the other case, only an insignificant reduction of the computation time or prohibitively large approximation errors In summary, efforts considering not only order reduction techniques but also efficient model formulations, cf.Popp et al. (2016bPopp et al. ( , 2017b)), as well as specialized numerical solvers, cf.Korolova et al. (2017), are important factors for a truly efficient simulation of nonlinear large scale interconnected dynamical systems.
Data availability.No data sets were used in this article.
Competing interests.The authors declare that they have no conflict of interest.Special issue statement.This article is part of the special issue "Kleinheubacher Berichte 2017".It is a result of the Kleinheubacher Tagung 2017, Miltenberg, Germany, 25-27 September 2017.

Figure 1 .
Figure 1.Stages of the hierarchical model formulation process within the mCCM framework; (a) compound component S containing individual components i ; (b) connected component C; (c) hierarchically connected component C.

Figure 2 .
Figure 2. Schematic representation of the modular concept of CoSimMA.

Figure 3 .
Figure 3. Example system for the mCCM modeling and TPWL reduction; (a) nonlinear chain network; (b) input section of the example circuit represented as mCCM component; (c) chain section of the example circuit represented as connected component; (d) interconnection scheme of the input section and the first two chain sections as mCCM representation of the considered example circuit.

Figure 4 .
Figure 4. Behavior of the weight functions w 1 . ..w 5 without units (p.u.) during the first second of simulation time.

Figure 5 .
Figure 5.Comparison of the original model with two reduced variants for the case 1.

Figure 6 .
Figure 6.Comparison of the original model with two reduced variants for the case 2.

Figure 7 .
Figure 7.Comparison of the Hankel singular values σ i for case 1 and case 2.