Basic principles and aims of model order reduction in compliant mechanisms

Abstract. Model order reduction appears to be beneficial for the synthesis and simulation of compliant mechanisms due to computational costs. Model order reduction is an established method in many technical fields for the approximation of large-scale linear time-invariant dynamical systems described by ordinary differential equations. Based on system theory, underlying representations of the dynamical system are introduced from which the general reduced order model is derived by projection. During the last years, numerous new procedures were published and investigated appropriate to simulation, optimization and control. Singular value decomposition, condensation-based and Krylov subspace methods representing three order reduction methods are reviewed and their advantages and disadvantages are outlined in this paper. The convenience of applying model order reduction in compliant mechanisms is quoted. Moreover, the requested attributes for order reduction as a future research direction meeting the characteristics of compliant mechanisms are commented.


Introduction
A new approach to develop a feed unit of small machine tools for small workpieces is based on the application of compliant mechanisms.Currently, non-intuitive design and optimization techniques are in progress as well as controlling, measuring and calibration strategies.To describe the mechanical behavior of a feed unit in an accurate way, very large and sparse finite element models arise.This leads to numerical simulations which require an unacceptable amount of time and memory space and motivates the introduction and application of model order reduction (MOR) methods.In this paper, MOR approaches are reviewed taking a first step towards its implemention in compliant mechanisms.To the authors' knowledge, applications and investigations of MOR in CMs are not existent.
In the present case, the feed unit consists of two major parts: (a) the compliant mechanism (CM) is a mechanical structure consisting of flexure hinges and rigid regions, including piezo-electric actuators and (b) the measure and control system supplying appropriate input signals for the mechanical structure via the actuators.Figure 1 illustrates a prototype of a CM appropriate to novel machine tools (the Correspondence to: M. Rösner (malte.roesner@hsu-hh.de)actuators are not shown), as described by Wulfsberg et al. (2010), embedded in the mechatronic system.
The utilization of CMs assuring the trajectory of the tool is constituted by their performance.CMs distinguish from previous standard machine tools being potentially more accurate, better scalable, cleaner, low-maintenance, less noisy and cheaper in manufacturing which makes them particularly suitable for small-scale applications.

System theoretical background
Assuming a linear time invariant system (LTI), the performance of the compliant mechanism can be specified by a system of second order ordinary differential equations (ODE) derived from linear finite element discretization, see e.g.Koutsovasilis and Beitelschmidt (2008).This system describes the reaction of the mechanical structure in answer to inputs from the control system and can be written in the time domain by means of the state space form in vectorial representation as: M e q(t) + D e q(t) + K e q(t) = B e u(t), y(t) = C e q(t), (1) where M e ∈ R n×n mass matrix, D e ∈ R n×n damping matrix, K e ∈ R n×n stiffness matrix, B e ∈ R n×p input matrix with u(t) ∈ R p input vector forming an external force vector, C e ∈ R m×n Published by Copernicus Publications.output matrix selecting the y(t) ∈ R m output vector of interest of the q(t) ∈ R n internal state vector corresponding to the structural characteristics of the CM.Moreover, n is the order of the system, also referred to as state space dimension.In case m = p = 1 and therefore B e and C e switch over to vectors as well as u(t) and y(t) to scalars, the systems is called single input -single output (SISO) otherwise multiple input -multiple output (MIMO).
To accurately reproduce the structural behavior, in most instances, the set of ODEs stated in Eq. ( 1) reaches a numerousness of degrees of freedom (DOF), stated by n in the present case.These large scale systems entail prohibitive computational time and memory space for calculation making them impractical for simulation, optimization, analysis and control.Therefore, MOR techniques are implemented to gain a low dimensional approximation of the high order system.MOR can be arranged in the simulation process, as shown in Fig. 2.
Depending on the reduction method, different representations of the system are required, which are illustrated below.
Several approximation methods rely on first order systems.The transformation of Eq. ( 1) to a first order descriptor system is realized by applying the state vector x(t) = [q T | qT ] T leading to the state-space dynamical system where A ∈ R n s ×n s system matrix, B ∈ R n s ×p input matrix, C ∈ R m×n s output matrix, E ∈ R n s ×n s descriptor matrix and n s = 2n.Since some reduction methods are not based on the representations written in Eqs. ( 1) and ( 2), respectively, but on the transfer function matrix H(s).Its derivation is described in the following.Using the Laplace transformation, the system specified in Eq. ( 1) is converted from time to frequency domain resulting in an algebraic system of equations given by where s is the complex-valued scalar parameter and Q(s), U(s), Y(s) are the Laplace transformed of q(t), u(t), y(t) and provided that homogeneous initial conditions q(t = 0) = q(t = 0) = 0 are existent.The transfer function matrix H(s) of the system can be specified by combining the two equations in Eq. ( 3), as

U(s),
(4) alternatively for the descriptor system in Eq. ( 2) where the transfer function matrix H(s) has a dimension of (m × p) in the MIMO case and, accordingly, is a scalar in the SISO case.The transfer function H(s) describes the correlation between input U(s) and output Y(s), disregarding the internal states of the system.Due to the fact that most MOR procedures are accomplished by means of projection, the approach is explained by a system represented in Eq. ( 2).The objective is to approximate the state vector x(t) ∈ R n s to a low-dimensional subspace.This is achieved by the substitution where V ∈ R n s ×q is the transformation matrix, x r (t) ∈ R q the reduced state vector, (t) the residual and q n s .Inserting the projection Eq. ( 6) into the linear system of ODEs in Eq. ( 2) results in a low-dimensional approximation EV ẋr (t) = AVx r (t) + BVu(t) + (t), y(t) = CVx r (t).( 7) Advantages Disadvantages -stability is preserved -not automatable -invariant to system representation -unnecessarily large reduced order models -physical information is preserved -not expandable to nonlinear systems -independent from the number of the systems inputs and outputs -no presetting of a minimal dimension of the reduced system possible -error estimation possible -computationally intense -implemented in software (Matlab, ANSYS) -inappropriate for large scale systems In general, this system is overdetermined having q unknowns but n s equations.To solve this problem and find a unique solution, Eq. ( 7) is pre-multiplied by a second transformation matrix W ∈ R n s ×q such that W T (t) = 0, the residual is equal to zero and The system in Eq. ( 8) is called general reduced order model (ROM) by projection.The number of inputs u(t) ∈ R p and outputs y(t) ∈ R m remains the same though the order of the system, expressed by the dimension of the reduced system matrices E r ∈ R q×q ,A r ∈ R q×q ,B r ∈ R q×p , C r ∈ R q×m and the reduced state vector x r (t) ∈ R q decreases from n s to q.
For most of the methods, the aim of model order reduction is to provide the projection matrices W and V such that their calculation is computationally efficient as well as automatable and the systems characteristics are preserved.The existence of a predefined error bound is desired in addition to the applicability in large-scale systems with an order up to a few hundred thousand.Details are specified e.g. in Eid (2009) and Rudnyi and Korvink (2006).

Reduction methods based on condensation
The reduction of a given system using condensation methods is the most commonly applied reduction method.Applying the time independent transformation matrices W T and V to the first equation in Eq. ( 1) gives the reduced system where the transformation matrices have to ensure the conservation of the potential and kinetic energy of the system.The reduced damping matrix D r may be modelled as linear combination of the mass and stiffness matrices (Rayleigh damping) with D r = αM r + βK r .
The construction of the transformation matrix V is carried out either by static, modal or mixed condensation which are explained briefly below, see Siedl (2008), Gasch and Knothe (1989), Bennini (2005) and Gugel (2009).In Table 1, the main advantages and disadvantages are listed for MOR using condensation-based reduction methods.

Static condensation
Static condensation is based on the partition of the systems degrees of freedom into dependent (MDOF -master degree of freedom, index m) and independent (SDOF -slave degree of freedom, index s) ones in a physical expedient way.This procedure is also the basis for sub-structure modeling and the generation of super elements.For the static case, q(t) = q(t) = 0, an external force vector is introduced by B e u(t) = F(t).The aforementioned partition into MDOF and SDOF applied to Eq. ( 1) results in The second row in Eq. ( 10) provides the so called static correlation between the MDOF and SDOF given by One obtains the reduced system by substituting q s (t) in the first row of Eq. ( 10) by Eq. ( 11).Rearranging ensues where Therewith, the transformation matrices in Eq. ( 8) have the following structure where I is the identity matrix.The identification of K red in this way is exact.In case of approximately reducing the mass matrix M e and damping matrix D e also with the transformation matrices in Eq. ( 13) from the static approach, the procedure is called Guyan reduction, as stated by Stelzmann et al. (2002) and Lienemann (2006).

Modal condensation
Using modal condensation, only a selected number of eigenmodes is taken into account.
where Φ ∈ R l×l is getting downsized by choosing a minor frequency range resulting in the reduced modal matrix Φ r ∈ R l×k .Using an orthogonal projection W = V as mentioned by Rickelt-Rolf ( 2009), the transformation matrices are chosen such that

Mixed condensation
Mixed condensation, known as reduction by Craig and Bampton (1968), combines static and modal reduction methods.In a first step, the system in Eq. ( 9) is similarly re-sorted with respect to SDOF and MDOF as described before.Secondly, the MDOF are blocked out by setting q m = 0 resulting in an auxiliary system of which the eigenmodes ϕ i are calculated subsequently.Analogue to the modal condensation, the reduced modal matrix Φ r is established by selecting certain eigenmodes as described by Siedl (2008).The SDOF is again linked to the MDOF via the static correlation and the transformation matrices are complemented to with Φ s = −K −1 ss K sm and Φ r = ϕ 1 , ϕ 2 , ••• , ϕ k as mentioned in Sect.3.1 and Sect.3.2, respectively.

Krylov subspace methods
Originally, Krylov subspace methods were developed for iteratively solving large and sparse linear systems of equations.
Multiple methods are available,

Calculate projection matrices
Choose numerical algorithm providing projection matrices Select appropriate Krylov subspace(s) Apply projection to the original system Original system Reduced system for example Generalized-Minimal-Residual-Method (GM-RES), Biconjugate-Gradient (BiCG,) Conjugate-Gradien-Square (CGS) and Transpose-Free-Quasi-Minimal-Residual (TFQMR), which can be found e.g. in Meister (2005), Kanzow (2005), Van der Vorst (2003) and Bai (2000) to name only a few.In case of MOR, Krylov subspaces are used to form the projection matrices V and W in Eq. ( 8).The general approach is based on the approximation of the transfer behavior of the original system by means of the transfer function H(s) in the frequency domain given in Eq. ( 4) and Eq. ( 5), respectively.MOR using Krylov subspaces is also well-known as moment matching, meaning that a certain number of moments of the reduced and the original system is equal.Hence, the moments of the transfer function and their matching in the sense of an approximation play a decisive role and are described below.Figure 3 illustrates the key steps using Krylov subspaces.
In Table 2, the main advantages and disadvantages are listed for MOR using Krylov subspaces, given in Salimbahrami (2006) and Eid (2009).

Transfer function and its moments
Considering the transfer function in Eq. ( 5) for the MIMO case, rearranging results in Taking into account the Neumann series, as described by Werner (2009) and Eid (2009), one obtains the following Taylor series This power series is called MacLaurin series, with M 0 j = C(A −1 E) j A −1 B called moments of the transfer function and its expansion point being s = 0. Depending on the expansion point, the moments and its corresponding MOR schemes are defined as: -Expansion point s = 0 results in a MacLaurin series with its moments defined by The corresponding MOR procedure is called Padé-Approximation.
-Expansion point s = s k results in a Taylor series with its moments defined by The corresponding MOR procedure is called Shifted Padé-Approximation, Rational Interpolation or Multipoint Padé-Approximation.
-Expansion point s → ∞ results in a Markov series with its moments defined by The corresponding MOR procedure is called Partial Realisation.
The users choice of the expansion point concerning its location and number in case of the Shifted Padé-Approximation, which is also suitable for multiple expansion points, affects the approximating characteristics.In the context of MOR, Padé-Approximation matches the behavior in the low frequency range whereas Partial Realisation fullfils in the high frequency range and Rational Interpolation in the user specified one.The usage of moment matching in sparse, large scale systems is due to the simplicity of the involved operations to compute these moments, namely matrix vector multiplication and matrix inversion.Details may be found in Eid (2009), Lehner and Eberhard (2006), Grimme (1997) and Salimbahrami (2006).

Krylov subspaces and moment matching
To provide the projection matrices V and W, Krylov subspaces are employed to match a certain number of the previously defined moments.A Krylov subspace is spanned by a succession of vectors as where A ∈ R n×n is a matrix, b ∈ R n is called starting vector and b,Ab,A 2 b,...,A q−1 b are the generated basic vectors spanning the q-th Krylov subspace.The first linearly independent basic vectors form the basis of this Krylov subspace.In case multiple starting vectors have to be considered, a q-th order block Krylov subspace can be rendered precisely by where A ∈ R n×n , B = [b 1 ,..., b p ] ∈ R n×p is the starting matrix containing the p linearly independent starting vectors.If p = 1, one obtains the standard Krylov subspace in Eq. ( 24) with the starting vector b.
For the reduced order system in Eq. ( 8), the basic vectors of a suitable Krylov subspace can be utilized to find the projection matrices V and W.An explicit calculation of the moments to match is not neccessary.As carried out in Koutsovasilis (2009) and Salimbahrami (2006) an input block Krylov subspace K q1 (A −1 E,A −1 B) and output block Krylov subspace K q2 (A −T E T ,A −T C T ) are to be used for an expansion point s = 0 to match some moments in Eq. ( 21). For If the colmuns of the matrix V form a basis of K q1 and K q3 , respectively, and the matrix W is chosen such that A r = W T AV is non-singular, then the first q1 m and accordingly q3 m moments match.A typical choice is W = V, named Galerkin (orthogonal) projection, otherwise W V Petrov-Galerkin (oblique) projection.Using only the input Krylov subspace is known as one-sided Krylov subspace method.Two-sided methods applying input and output Krylov subspaces match q1 m + q2 m and accordingly q3 m + q4 m moments.To calculate the desired projection matrices V and W, a wide range of numerical, mostly iterative algorithms is avaiable.The most popular ones are Arnoldi algorithm, twosided Arnoldi algorithm and Lanczos algorithm which are not further investigated in this paper.Details and further algorithms can be found in Bechtold et al. (2007), Meister (2005), Grimme (1997), Krohne (2007), Salimbahrami (2006), Antoulas et al. (2001), Eid (2009) and Bai (2002).

Control theory methods based on singular value decomposition
Methods using the singular value decomposition (SVD) aim for an approximation of matrices by lower ranked matrices, as carried out by Hackbusch (2009).The reduced order model (ROM) is found by deleting those system states that are both least controllable and observable.Figure 4 gives an overview of the key steps using singular value decomposition.In Table 3, the main advantages and disadvantages are listed for MOR using SVD methods.

Controllability and observability of a system
Considering the system in Eq. ( 2) represented by the two important control theory quantaties controllability P and observability Q can be introduced, as done by Lienemann (2006), via the Gramians and The concept of controllability examines how the states x(t) of a system are linked to its inputs u(t) whereas observability deals with the connection between the states x(t) and its outputs y(t), as delineated by Eid (2009).Under certain conditions, described e.g. by Bechtold et al. (2007) and Koutsovasilis (2009), the Gramians can be found by solving the two system related Lyapunov equations given by

The singular value decomposition
As known from linear algebra and described by Gugel (2009) and Hackbusch (2009), any matrix A ∈ R m×n can be decomposed into the product of three matrices, more precisely an orthogonal matrix U, a diagonal matrix Σ and the transpose of an orthogonal matrix V such that where the matrix  30) is called a singular value decomposition (SVD) of A and provides a basis for finding a lower ranked matrix A r which best approximates A such that The task is to find an optimal approximation in a certain system norm, eg. the spectral norm (H 2 -norm), as specified by Antoulas and Sorensen (2001), Hackbusch (2009) and Dahmen and Reusken (2008), via the Schmidt-Mirsky theorem which results in the minimisation of the approximation error such that min rank(A r )≤rank(A) In the following, the most commonly used SVD procedures are specified, namely balanced truncation approximation, singular perturbation approximation and Hankel norm approximation.They distinguish in the rating of the approximation error and the procedure of its minimization, see Benner et al. (2004).

Balanced truncation approximation
In case of a stable system, balanced truncation approximation (BTA) assumes that the system is balanced first and afterwards truncated.Based on the system representation given in Eq. ( 26), the balanced partition of the system matrices and state vector can be quoted as with Â11 ∈ R r×r , B1 ∈ R r×m and Ĉ1 ∈ R p×r .At all times, a balancing transformation is feasible for minimal order systems.The McMillan degree of the system is specified as the order of any minimal state-space realization of the transferfunction matrix, as stated by Sasane (2002) and Baur and Benner (2008).
The singular values are the eigenvalues of the observability P and controllability Q Gramians: σ i = λ i (PQ) which are called Hankel singular values (HSV) of Σ, e.g.described by Gugercin and Antoulas (2004) and Antoulas and Sorensen (2001).The Gramians are equal and diagonal.The balanced truncation approximation is predicted on transforming the state-space-system into a balanced realization.Those state variables that are least observable and controllable, and therefore related to the smallest Hankel singular values, are truncated with respect to the error bound between original and reduced transfer function The ROM obtained by BTA is related to as mentioned by Obinata and Anderson (2000).BTA involves an approximation error in the low-frequency region as mentioned by Bechtold et al. (2007), Obinata and Anderson (2000) and Baur and Benner (2008).For details and aspects in numerical calculation see e.g.Benner and Quintanaort'i (2004), Bechtold et al. (2007), Koutsovasilis (2009), Antoulas (2005) and Gugercin and Antoulas (2004).

Singular perturbation approximation
The singular perturbation approximation (SPA) is closely related to the balanced truncation approximation and ensures zero error at zero frequency and is applicable to nonlinear systems.Refering to the expression in Eq. ( 34) and as explained by Liu and Anderson (1989) and Antoulas et al. (2001) the reduced system is identified by with the reduced system matrices Ã . In contrast to BTA, SPA has an approximation error at high frequencies, but it matches at low frequencies, see Obinata and Anderson (2000).

Hankel norm approximation
The Hankel norm approximation (HNA) employs the socalled Hankel norm which is defined as the maximal Hankel singular value of the system given in Eq. ( 26): The aim is to minimize the approximation error between the original system transfer function H(s) and the reduced one H r (s) with order r via the Hankel norm as stated by Glover (1984), Green and Limebeer (1995), Bechtold et al. (2007), Benner et al. (2004), Antoulas (2005) and Antoulas and Sorensen (2001) where also the computation of the optimal Hankel norm approximation is performed.In Eq. ( 39) Hr (s) denotes all transfer functions of McMillan degree less than or equal to r.

Conclusion: benefits and aims for model order reduction in compliant mechanisms
In this work, system theoretical basics are reviewed which constitute the origin for model order reduction (MOR).Three different procedures, namely modal reduction, Krylov subspaces and singular value decomposition (SVD) based methods, are described with their associated characteristics.The application of MOR can be motivated by challenges arising in the numerical characterisation of compliant mechanisms (CM) which can be traced back to the numerousness of ODEs required to specify their behavior in an accurate way.In most cases, large finite element models arise to capture a detailed representation of the geometrical domain and the dynamical performance of the CM.This circumstance is unfeasible, particularly with regard to further investigations such as simulation, optimization and control.For this reason, MOR in CMs seems to be a purposeful approach to generate efficient models and would benefit the design and development as well as analysis process of CMs.
The specific feature of CMs exists in their mechanical topology as they consist of rigid parts connected by flexure hinges.This characteristic has effects on the potentiality of MOR and requires further investigations to provide a suitable algorithm.Especially, the algorithm has to notably incorporate the distinct differences between the compliant and rigid parts of the structure with regard to their mechanical properties and influence on the systems transfer behavior.Moreover, having a small error between original and reduced system implying a global error bound, the procedure should preserve the passivity as well as stability, has to be automatable, stable and computational efficient also for high order systems.
In future work, different MOR schemes applied to CM will be investigated and reasonable combinations among them will be analysed to benefit from their advantages.Furthermore, a specific adjustment and upgrading is desired.Out of this, the aim is to provide a reduction procedure meeting the special circumstances of CMs.

Figure 3 .
Figure 3. Key steps of dimensional reduction using Krylov subspaces.

Figure 4 .
Figure 4. Key steps of dimensional reduction using singular value decomposition (SVD).
values σ i = λ i (A T A) ≥ σ i+1 ≥ 0 on its diagonal in descending order with rank(A) =k and k = min{m,n}.The left singular vectors U = [u 1 u 2 ... u m ] ∈ R m×m are the orthonormal vectors of A A T , the right singular vectors V = [v 1 v 2 ... v n ] ∈ R n×n are the orthonormal vectors of A T A. The term in Eq. (

Table 1 .
Characteristics of condensation-based methods.

Table 2 .
Characteristics of Krylov subspace methods.

Table 3 .
Characteristics of SVD methods.