Fundamentals of Model Order Reduction

Massachusetts Institute of Technology Department of Electrical Engineering and Computer Science 6.245: MULTIVARIABLE CONTROL SYSTEMS by A. Megretski Fundamentals of Model Order Reduction1 This lecture introduces basic principles of model order reduction for LTI systems, which is about finding good low order approximations of high order systems. 8.1 Setting Up A Model Order Reduction Problem This section describes a number of different ways of setting up a model reduction problem. In contrast with the case of LTI feedback optimization, when the most natural problem formulation has a satisfactory solution (H2 or H-Infinity optimization via Riccati equa­ tions), no efficient solution is known for the most natural formulation of the model order reduction problem. This leads to a variety of compromise approaches. 8.1.1 An abstract setup To define a dynamical system model complexity reduction problem, one typically needs (at least) the following two ingredients. (a) A numerical measure of complexity, applicable to system models under considera­ tion. For example, complexity of an LTI system can be measured by its order, while complexity of a finite state automata can be measured by the number of its states. (b) A numerical measure of distance between two systems. For example, the L2 norm of the difference of impulse responses can serve as a measure of distance between two LTI systems. 1 Version of March 8, 2004 2 A typical model reduction problem will then be formulated as the task of finding a ”re­ duced” system Ŝ = Ŝk of complexity not larger than a given threshold k, such that the distance between Ŝ and a given ”complex” system S is as small as possible. Alterna­ tively, a maximal admissible distance between S and Ŝ can be given, in which case the complexity k of Ŝ is to be minimized. As is suggested by the experience of robust control, useful measures of distance between S and Ŝ can be defined as the induced gain (for example, the L2 gain) from a testing input f to the output matching error e, as shown on Figure 8.1. In this case, for design f S + � �− e Ŝ Figure 8.1: Comparing S and Ŝ or analysis purposes, S can be represented as a series connection (i.e. a sum) of Ŝ and an ”uncertain” error system � = S − Ŝ. 8.1.2 H-Infinity Optimal Model Reduction Setup For LTI systems G given in an input/output format (for example, by their transfer matri­ ces G = G(s) or by the impulse response g = g(t)), complexity is well represented by the system order, defined as the number of states in a state-space realization which is both controllable and observable (also called the McMillan degree). The H-Infinity model reduction problem is that of finding a stable LTI transfer matrix ˆ ˆ k of order smaller than a given number k such that ≡W −1 (G − G ˆ k )≡� is as small G=G as possible, where G, W are given stable transfer matrices (W −1 is also assumed to be stable), and ≡�≡� denotes H-Infinity norm (L2 gain) of a stable system �. As a result of model order reduction, G can be represented as a series connection of a lower order “nominal plant” Ĝ and a bounded uncertainty (see Figure 8.2). In most cases, (even when W ≥ 1) the H-Infinity optimal model order reduction is a problem without a known good solution. This can be explained (informally) as a consequence of non-convexity of the set of all stable transfer matrices of order less than k. In particular, one can expect many local minima in an attempt at a gradient-based optimization of ≡W −1 (G − Ĝ)≡� . 3 G(s) � � + ˆ G �(s) W (s) Figure 8.2: Model order reduction and uncertain models 8.1.3 Approximation by Linear Combinations vs. Model Reduction A much easier problem is that of finding the best approximation of a given LTI system ˆ of given LTI systems Gk : G by a linear combination G ⎦ ˆ � min, G ˆ= ≡G − G≡ c r Gr , where ≡ · ≡ is some system norm. For many approximation quality measures, this problem has a relatively simple solution. For example, if ≡ · ≡ = ≡ · ≡2 is the H2 norm then the optimization becomes a standard least squares problem reducible to solving a system of linear equations. If ≡ · ≡ = ≡ · ≡� is the H-infinity norm, the optimization is reduced to solving a system of Linear Matrix Inequalities (LMI), a special class of convex optimization problems solvable by an efficient algorithm, to be discussed later in the course. While the technique of approximation by linear combinations should not be expected to yield a close-to-optimal relation between the order of Ĝ and the size of model matching error ≡G − Ĝ≡, it is frequently used in the initial stage of a model reduction process as a tool for approximating infinite dimentional systems by finite order models. In particular, a least squares algorithm can be used to derive a high quality approximation G0 of a nonrational transfer matrix G, such that G0 has a high but finite order. Then, if ≡G − G0 ≡� is small enough, model reduction of G can be replaced by model reduction of G0 . 8.1.4 Hankel Model Reduction The so-called Hankel optimal model reduction algorithm is based on finding a stable m-th ˆ of a given stable LTI system G which minimizes the so­ order LTI approximation G 4 called Hankel norm ≡G − Ĝ≡H of model reduction error. Technical details of definitions and optimization algorithms associated with the Hankel norm will be discussed in the following sections. For the user, it is important to know that Hankel norm ≡�≡ H of a stable system � is never larger than its H-Infinity norm ≡�≡� , hence solving the Hankel norm optimal model reduction problem yields a lower bound for the minimum in the HInfinity norm optimal model reduction. Moreover, H-Infinity norm of model reduction error associated with a Hankel norm optimal reduced model is typically close to this lower bound. Thus, Hankel norm optimal reduced model can serve well as an H-Infinity suboptimal reduced model. Hankel optimal model reduction appears to be one of the best model reduction options available. However, this approach also has some weak spots. One problem is the cubic (with respect to the order of the original system G) growth of complexity of the model reduction algorithm, which makes it impractical to apply it to systems of order 10 3 and larger. Another (perhaps related) weakness lies in the difficulty of adapting the algorithm to a situation when G is not perfectly known, and, alternatively, samples of its time or frequency response are available. 8.1.5 Non-Optimizing Model Reduction Due to the prohibitively high cost of Hankel model reduction, applications dealing with very large order systems frequently resort to numerically inexpensive methods not promis­ ing any optimality at all. We will refer to them as non-optimizing model reduction tools, stressing that no model reduction error is minimized when the methods are applied. In this subsection we introduce some of these methods briefly, using the SISO model reduc­ tion, in which a given stable transfer function G = G(s) is to be approximated by the ˆ ratio G(s) = p(s)/q(s), where p, q are polynomials of order m. One popular approach is moments matching. In the simplest form of moments matchˆ ing, an m-th order approximation G(s) = p(s)/q(s) (where p, q are polynomials of order m) of a SISO transfer function G(s) is defined by matching the first 2m + 1 coefficients of a Taylor expansion 1 G(s) = G(0) + G� (0)s + G�� (0)s2 + . . . 2 ˆ This way, one gets a system of 2m + 1 by the coefficients of a similar expansion of G. equations with 2m + 1 parameters. There is a whole branch of model reduction science of writing these equations in a numerically robust form, of using other (or multiple) frequencies for matching Taylor expansion coefficients, etc. However, the following major flaw remains with the approach: it does not guarantee any degree of near-optimality of the 5 ˆ In practice, the moments matching reduced models, nor even stability of the resulting G. approach looses dramatically to Hankel and other optimization-based model reduction algorithms. Another popular method of model reduction is based on minimizing the error of match­ ing a fixed output. For example, assume that samples yk = y(kT ), k = 1, 2, . . . at a fixed sampling rate T > 0 are available for the impulse response of a given CT LTI system. If the system has order m, there would exist a Schur polynomial q(z) = qm z m + qm−1 z m−1 + · · · + q1 z + q0 , with qm ∞= 0, such that qm yk+m + qm−1 yk+m−1 + · · · + q1 yk+1 + q0 yk = 0 ˆ for all k > 0. The idea is to define the denominator of the reduced model G(s) in terms of a polynomial q = q(z) which minimizes the sum of squares of ek = qm yk+m + qm−1 yk+m−1 + · · · + q1 yk+1 + q0 yk subject to a normalization constraint imposed on the coefficients of q (otherwise the minimum is achieved at q ≥ 0). When the normalization is given by qm = 1, we get a classical least squares system identification algorithm. When the normalization is given by m ⎦ |qk |2 = 1, k=0 we get a version of the proper orthogonal decomposition algorithm. The trouble with these approaches is very much the same as with moments matching: no reason to ex­ pect near-optimality of the reduced model, and, unless significant performance-reducing modifications are introduced, the reduced model generated by these algorithms is not necessarily stable. 8.2 Hankel Operators and System Balancing Here we introduce the fundamental notions of Hankel operator and Hankel singular num­ bers associated with a given stable LTI system. For practical calculations, balancing of stable LTI system is defined and explained. 6 8.2.1 Hankel Operator The “convolution operator” f ∈� y = g � f associated with a LTI system with impulse response g = g(t) has infinite rank whenever g is not identically equal to zero. However, with every LTI system of finite order, it is possible to associate a meaningful and repre­ sentative linear transformation of finite rank. This transformation is called the Hankel operator. A Hankel operator HG can be associated with every LTI system G which has a finite L2 gain. It is important that G does not have to be causal. The domain D(HG ) of HG consists of all fignals f = f (t) of finite energy such that f (t) = 0 for t → 0. The result of applying the Hankel operator HG to a signal f ≤ D(HG ) is a signal y = y(t) which, for t > 0, equals the response of G to input f , and equals zero otherwise. In other terms, HG : L2 (−∗, 0) ∈� L2 (0, ∗) is defined by ⎡ � (HG f )(t) = g(t + ψ )f (−ψ )dψ, 0 where g = g(t) is the impulse response of G. Let G = G(s) be a m-by-k matrix-valued function of complex argument which is defined on the extended imaginary axis jR � {∗}, where it satisfies the conditions of real symmetry (i.e. elements of G(−jρ) are the complex conjugates of the corresponding entries of G(jρ)) and continuity (i.e. G(jρ) converges to G(jρ0 ) as ρ � ρ0 for all ρ0 ≤ R and for ρ0 = ∗). Note that a rational function G = G(s) with real coefficients satisfies this condition if and only if it is proper and has no poles on the imaginary axis. Let Lk2 (−∗, 0) denote the set of square integrable functions f : R ∈� Rk such that m f (t) = 0 for t → 0. Let Lm 2 (0, ∗) be the set of square integrable functions g : R ∈� R such that g(t) = 0 for t ∪ 0. The Hankel operator associated with G is a map H G : Lk2 (−∗, 0) ∈� Lm 2 (0, ∗) defined by ⎤ g0 (t), t > 0, (HG f )(t) = g(t) = 0, otherwise, where g0 = g0 (t) is the inverse Fourier transform of g̃0 (jρ) = G(jρ)f˜(jρ). Less formally, the Hankel operator HG associated with an LTI system G takes an input f which equals zero for for t → 0, and produces a signal g which equals the corresponding output of G for t > 0, and equals zero for t ∪ 0. 8.2.2 Rank and Gain of a Hankel Operator As a linear transformation HG : Lk2 (−∗, 0) ∈� Lm 2 (0, ∗), a Hankel operator has its rank and L2 gain readily defined: the rank of HG is the maximal number of linearly independent 7 outputs (could be plus infinity), the L2 gain is the minimal upper bound for the L2 norm of HG f subject to the L2 norm of f being fixed at one. Remember that the order of a rational transfer matrix G is defined as the order of its controllable and observable realization. The order of a non-rational transfer matrix equals infinity. The L-Infinity norm of a transfer matrix is defined similarly to the H-Infinity norm, but without the requirement of stability: ≡G≡� = sup λmax (G(jρ)). ��R There is a simple but very useful relation between the rank and norm of a Hankel operator HG , on one hand, and the order and H-Infinity norm of the original LTI system G, on the other. Theorem 8.1 For all transfer matrices G, the rank of HG does not exceed the order of G, and the L2 gain of G does not exceed the L-Infinity norm of G. Moreover, for stable transfer matrices G, the rank of HG equals the order of G. In other words, rank of HG equals the order of the stable part of G, and L2 gain of HG is never larger than ≡G≡� . Proof To show the gain/L-Infinity norm relation, let us return to the definition of H G in the previous subsection. Note that the L2 norm of g0 is not larger than ≡G≡� ≡f ≡2 . On the other hand, ≡g≡2 ∪ ≡g≡2 . To show that the rank of HG equals the order of the stable part of G, note first that the unstable part of G does not have any effect of HG at all. Then, for a stable G, the map of f into the inverse Fourier transform g0 of g̃0 = Gf˜ is a causal linear system. Hence, the only way for a past signal f ≤ L2 (−∗, 0) to affect the future of the output is through the system’s state at zero. Naturally, the so-called Hankel norm ≡G≡H of G is defined as the L2 gain of the corresponding Hankel operator HG . 8.2.3 Matrix rank reduction The common feature matrices and linear systems share is that both describe linear trans­ formations. Before addressing the problem of LTI system model order reduction, it is instructive to consider the case of matrices, in which a similar task is defined as the matrix rank reduction problem: given a matrix M , and an integer k > 0 find a matrix ˆ =M ˆ k of rank less than k which minimizes λmax (M − M ˆ ). M 8 Since the set of matrices of given dimensions of rank less than k is not convex (unless k is larger than one of the dimensions), one can expect that the matrix rank reduction prob­ lem will be difficult to solve. However, in fact it has a relatively simple computationally efficient solution. Let m ⎦ M= ur λr vr� r=1 be a singular value decomposition of M , which means that the families {ur }m r=1 and m m {vr }r=1 are orthonormal, and {λr }r=1 is a monotonically decreasing sequence of positive numbers, i.e. ⎤ 1, i = r, � � λ1 → λ2 → · · · → λm > 0. ui ur = v i v r = 0, i ∞= r, Here m is the rank of M , the numbers λr are called singular numbers of M , and the vectors vr , ur are, respectively, right and left singular vectors of M . The vector v1 has unit length |v1 | = 1 (the Euclidean norm) and yields a maximal (over the set of all vectors of unit length) length λ1 = |M v1 | when transformed with M . Vector u1 is defined by M v1 = λ1 u1 . The vector v2 has unit length, is orthogonal to v1 , and yields a maximal (over all vectors of unit length and orthogonal to v1 ) length λ2 = |M v2 | when transformed with M . Vector u2 is defined by M v2 = λ2 u2 . In general, the vector vr has unit length, is orthogonal to v1 , . . . , vr−1 , and yields a maximal (over all vectors of unit length and orthogonal to v1 , . . . , vr−1 ) length λr = |M vr | when transformed with M . Vector ur is defined by M vr = λr ur . Another useful interpretation of singular vectors vr , ur and singular numbers λr is by the eigenvalue decompositions (M � M )vr = λr2 vr , (M M � )sur = λr2 ur . ˆ =M ˆ k of rank less than k which minimizes λmax (M − M ˆ ) is given An approximation M by ˆk = M k−1 ⎦ ur λr vr� r=1 and yields an approximation error of λmax (M − M̂k ) = λk . In other words, matching the original linear transformation M on the first k − 1 inputs of highest amplification yields ˆ. the best reduced rank approximation M 9 ˆ . For example, when In most cases there are many optimal approximations M ⎣ � 4 0 0 M =� 0 2 0 � 0 0 1 then both � ⎣ � 4 0 2 4 0 0 M̂a = � 0 0 0 � , M̂b = � 0 0 0 � 2 0 1 0 0 0 ⎣ are optimal rank one approximations of M . The framework of optimal matrix rank reduction can be easily extended to the class of linear transformations M from one Hilbert space to another. (In the case of a real nby-m matrix M , the Hilbert spaces are Rm and Rn .) One potential complication is that the vectors vr of maximal amplification do not always exist. In these notes, the following definition of singular numbers λr = λr (M ) will be used: if the first k right singular vectors v1 , . . . , vk of M can be defined, but the supremum � = sup{|M v| : |v| = 1, v orthogonal to v1 , . . . , vk } is not achievable, then λr (M ) = � for all r > k. 8.2.4 SVD of a Hankel operator Any causal stable LTI system with a state-space model dx/dt = Ax + Bf, y = Cx + Df defines a map past input ∈� x(0) To produce x(0) = x0 , the energy of the past input must be at least x�0 �0 , where �0 is an arbitrary solution of Wc �0 = x0 , ⎡ � � Wc = eAt BB � eA t dt > 0 0 is the controllability Grammian of the system. (If Wc is not positive definite, and equation Wc �0 = x0 is not solvable with respect to �0 , x0 is an unreachable state.) The lowest energy input f (t) producing x0 is given by � f (−t) = B � eA t �0 (t > 0). 10 Note that calculation of Wc is easy via the Lyapunov equation AWc + Wc A� = −BB � The energy of the future output produced by the initial state x(0) = x0 , provided zero input for t > 0 equals x�0 Wo x0 , where ⎡ � � W0 = eA t C � CeAt dt 0 is the observability Gamian of the system. The output produced by x0 is given by y(t) = CeAt x0 , and calculation of Wo is easy via the Lyapunov equation Wo A + A� Wo = −C � C SVD of a Hankel operator H can be expressed in terms of its Gramians: 1/2 1/2 Let wi be the normalized eigenvectors of R = Wc Wo Wc , i.e. Rwi = �i wi , �1 → �2 → . . . , �m > 0, �m+1 = 0 The SVD of H is given by H= m ⎦ uk λk vk� , k=1 1/2 where λk = �k , uk (t) = CeAt Wc1/2 wk �k −1/2 vk (t) = 8.2.5 (t > 0) � −A� t �−1 Wo Wc1/2 wk k B e (t < 0) Technical details of the proof Let M : Rn ∈� L2 (0, ∗) be defined by (M x0 )(t) = CeAt x0 and let N : L2 (−∗, 0) ∈� Rn be N u(·) = ⎡ 0 e−At Bu(t)dt −� 11 By the definition of Wc , Wo , Wo = M � M, Wc = N N � Hence M, N can be represented in the form M = U Wo1/2 , N = Wc1/2 V � where linear transformations U, V preserve the 2-norm. Since the Hankel operator under consideration has the form H = M N = U Wo1/2 Wc1/2 V � in order to find SVD of H, it is sufficient to find SVD of F = Wo1/2 Wc1/2 Since F � F = R, the SVD is given by Wo1/2 Wc1/2 = m ⎦ 1/2 u¯k �k v¯k� k=1 where −1/2 uk = Wo1/2 Wc1/2 wk �l , vk = Wc1/2 Wo p1/2 wk �−1 k = wk .