- Email: [email protected]
ℋ2 Model Reduction by Means of Gröbner Bases and its Extension to the Parametric Case
16th IFAC Symposium on System Identification The International Federation of Automatic Control Brussels, Belgium. July 11-13, 2012
H2 Model Reduction by Means of Gr¨ obner Bases and its Extension to the Parametric Case Masaaki Kanno ∗ ∗
Institute of Science and Technology, Academic Assembly, Niigata University, Niigata, 950-2181, Japan (e-mail: [email protected]).
Abstract: In this paper a solution approach for parametric H2 model reduction based on Gr¨obner bases is presented. The approach utilizes a na¨ıve parametrization of the approximant, which leads to a rather simple set of algebraic equations that is solved by means of Gr¨obner bases. This allows one to deal with the parametric case, and analysis of the optimal solution in the presence of parameters becomes possible. Numerical examples show that the optimal approximant may change discontinuously as a parameter value varies, also showing the existence of systems with two optimal approximants. Keywords: H2 model reduction; algebraic methods; parametric case; multiple optimal solutions 1. INTRODUCTION Modelling is instrumental in modern science and engineering. A good model elucidates the nature of the object under investigation and helps one to construct a system that interacts with it. The meaning of a ‘good’ model can vary. An accurate model that reproduces the behaviour of the true object precisely may be beneficial in simulation, but can conceal beneath too much detail the real working mechanisms that want examination (Anderson, 1978). In modern control system design techniques, a detailed and complicated model may result in a controller of high orders, which can cause implementation issues due to computational complexity and numerical difficulties and may also be fragile against perturbation. The problem of model reduction has long been a subject of research in the systems and control field (Obinata and Anderson, 2001). In order to derive a simplified model from a complicated one, numerous model reduction techniques have been developed. One of popular methods is balanced truncation initiated by Moore (1981), which is further advanced by frequency-weighting (Ghafoor and Sreeram, 2008). Another popular approach is the Hankel norm approximation method (Glover, 1984). The advantage of those approaches is that they admit constructive algorithms. By setting some cost functions to measure the ‘closeness’ of the original system and the simplified one, model reduction can be formulated as optimization problems. Solution approaches appear for H2 and H∞ model reduction, e.g., Skelton et al. (1998). Another optimization problem setting relevant in feedback control is ν-gap model reduction investigated in Cantoni (2001) since the ν-gap is an appropriate metric to measure the closeness of two systems in a feedback loop. This paper considers H2 model reduction. The use of the H2 -norm to evaluate the fidelity of the reduced order model is advantageous in that this criterion admits both the time-domain and frequency-domain interpretations. 978-3-902823-06-9/12/$20.00 © 2012 IFAC
710
Investigation on this problem has a long history, and a wide range of techniques have been proposed. Many approaches rely on the first order necessary condition for optimality (Meier III and Luenberger, 1967) and its extension to multivariate case (Gugercin et al., 2008). A recent result numerically finds approximants for systems with thousands of states effectively (Antoulas, 2005; Gugercin et al., 2008). The problem of H2 model reduction is a nonconvex problem and it is intrinsically difficult to find the global optimal approximant. Even though a numerically effective approach can deal with a very large-scale system, the quality of the resulting approximant thus depends on the choice of the initial guess. The approach employed in this paper is based on an algebraic approach, more specifically Gr¨obner bases (Cox et al., 2007). The advantage of such approach over typical numerical approaches is that it may achieve global optimality and also can handle the parametric case. The drawback (and a major one) is computational complexity due to the reliance of symbolic manipulation and exact arithmetic. The aims of the paper is first to present the current state of an algebraic approach and further to investigate the H2 model reduction problem with the aid of such approach suggesting how difficult the problem can be. The paper is organized as follows. Section 2 states the problem formulation, followed by the suggested solution approach in Section 3. Three numerical examples are given in Section 4, where the approach is demonstrated and further it is shown how the number of candidates and the optimal reduced order model may change as a parameter value changes. Also systems with two global optimal solutions are presented. Some concluding remarks are made in Section 5. 2. PROBLEM FORMULATION The problem of H2 model reduction considered in this paper is formulated as follows. Given an n-th order stable 10.3182/20120711-3-BE-2027.00160
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
single-input, single-output system G(s), the task is to find ˆ an m-th order stable system G(s) (where m < n) that minimizes
ˆ ,
G(s) − G(s) 2 i.e., to find the best approximant in the H2 sense. The problem has been studied for decades and a wide variety of mathematical tools have been applied to solve the problem. An algebraic approach is utilized in this paper and one of early work in that direction is Hanzon and Maciejowski (1996). 3. SOLUTION APPROACH The solution approach which is employed here and initiated in Kanno (2010) relies on an algebraic method. The coefficients of the denominator of the approximant transfer function and those of an auxiliary polynomial are taken as unknown variables, and a set of algebraic equations in terms of those variables are derived from an optimality condition. Then a set of solutions, which constitutes a set of candidates, are calculated by means of an important algebraic method called Gr¨obner bases (Cox et al., 2007). Of such candidates, a stable system with real coefficients that achieves the minimal H2 -norm is the true optimal approximant. The main difference of the approach employed here from that of Hanzon and Maciejowski (1996) is a simple parametrization of the approximant. The approach by Hanzon and Maciejowski (1996) exploits the Schwarz-like canonical form which simplifies the stability criterion. The approach in this paper requires an extra effort and depends upon, e.g., the Routh-Hurwitz criterion to guarantee the stability of the approximant. Nevertheless simpler equations are obtained and the approach can deal with the parametric case with relative ease, which allows one to appreciate the strength of the algebraic approach. For instance, one may analyse how the optimal approximant changes as parameter values vary. Write the original system G(s) and its reduced order ˆ approximant G(s) to be sought as ˆ (s) N N (s) ˆ , G(s) = , G(s) = ˆ D(s) D(s) ˆ where D(s) and D(s) are stable monic polynomials of ˆ (s) are degree n and m, respectively, and N (s) and N polynomials of degree strictly less than n and m, respecˆ tively. For the optimal approximant G(s), there exists an auxiliary polynomial Q(s) of degree at most n−m−1 that satisfies the following equation (Meier III and Luenberger, 1967): 2 ˆ ˆ ˆ (s)D(s) . Q(s) = N (1) N (s)D(s) − D(−s) In this paper it is a standing assumption that the optimal ˆ reduced order model can be obtained by finding D(s) and Q(s) satisfying (1). Write ˆ D(s) = sm + am−1 sm−1 + · · · + a1 s + a0 , Q(s) = bn−m−1 sn−m−1 + · · · + b1 s + b0 , and compute the remainder R(s) on dividing the left hand side of (1) by D(s). Since R(s) has to be identically 0, its
711
coefficients, which are polynomials in ai ’s and bj ’s, have to be all zero. One thus obtains a set of algebraic equations. ˆ (s) can be computed Once the equations are solved, N by (1). In order to solve the obtained equations, a Gr¨obner basis of the polynomial parts of the equations is computed, which yields an equivalent set of equations that is easier to ‘solve’. It is further assumed that the system of equations has a finite number of solutions, i.e., the solution set is 0-dimensional. The Gr¨obner basis with respect to the lexicographic ordering bn−m−1 ≻ bn−m−2 ≻ · · · ≻ b0 ≻ am−1 ≻ am−2 ≻ · · · ≻ a0 is generically in the shape form: { Sf (a0 ) , a1 − h1 (a0 ) , . . . , am−1 − hm−1 (a0 ) , b0 − hm (a0 ) , . . . , bn−m−1 − hn−1 (a0 ) } . The set of zeros is then obtained by firstly finding the real roots of Sf (a0 ) and then substituting the values in the expressions hi (a0 ). Some of the zeros may yield unstable systems, which are to be discarded at this point. One then ˆ computes the H2 -norm of the error system G(s) − G(s) for each candidate and picks up the candidate giving the smallest H2 error, which is the optimal approximant. In the parametric case, Sf (a0 ) contains parameters, namely, Sf (a0 ) is a polynomial in a0 whose coefficients are functions in parameters. If parameters enter in polynomial form in the coefficients of G(s), then Sf (a0 ) will also be a polynomial in parameters. Once a parametric Sf (a0 ) is obtained, one can exploit such expression to compute the optimal approximant for each set of parameter values, and thus can utilize the result of off-line computation and avoid recomputation from the start. Namely one has to substitute the parameter with the value and compute positive roots of a polynomial. Then a finite number of candidates are obtained, and the H2 -norm of the error system for each case is to be calculated for determining the true optimal approximant. Further analysis may provide information on how many candidates can be obtained. This point is demonstrated in numerical examples in the subsequent section. 4. NUMERICAL EXAMPLES 4.1 A 3rd Order System The following 3rd order system with a parameter ζ is approximated by a 1st order system: 1 50 G(s) = + 2 s+1 s + 20ζs + 100 s2 + (20ζ + 50)s + 150 N (s) = 3 =: , 2 s + (20ζ + 1)s + (20ζ + 100)s + 100 D(s) where ζ > 0 (stability requirement). The denominaˆ ˆ tor D(s) of the approximant G(s) and the auxiliary polynomial Q(s) are written as ˆ D(s) = s + a0 , Q(s) = b1 s + b0 . (2) The quotient and the remainder on dividing the left hand side of (1) by D(s) are
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
1 − b1 , 2a0 b1 + a0 − b0 + (20ζ + 1)b1 + 49 s2 + −a20 b1 + 2a0 b0 + (20ζ + 50)a0 + (20ζ + 100)b1 −20ζ + 50 s − a20 b0 + 150a0 + 100b1 − 100 , respectively. The remainder has to be identically 0, and thus what is to be found is the zeros of the following ideal: {2a0 b1 + a0 − b0 + (20ζ + 1)b1 + 49 , − a20 b1 + 2a0 b0 + (20ζ + 50)a0 + (20ζ + 100)b1 − 20ζ + 50 , − a20 b0 + 150a0 + 100b1 − 100} . (3) Computing the Gr¨obner basis with respect to the lexicographic ordering b1 ≻ b0 ≻ a0 , one can get the Gr¨obner basis in shape form for the generic case: {a50 + (40ζ + 149)a40 + (400ζ 2 + 960ζ + 500)a30 + (−400ζ 2 + 6000ζ − 5050)a20 − 3000ζa0 − 15000 , 200(ζ 2 − 1)(20ζ − 101)2 b1 + (−400ζ 2 + 20ζ + 299)a40 + (−16000ζ 3 − 58600ζ 2 + 15950ζ + 44151)a30 + (−160000ζ 4 − 368000ζ 3 + 8600ζ 2 + 433530ζ + 80200)a20 + (240000ζ 4 − 1824000ζ 3 + 2910600ζ 2 − 3135450 + 1821000ζ)a0 − 80000ζ 4 + 1408000ζ 3 + 979800ζ 2 − 1256500ζ − 929800 , 200(ζ 2 − 1)(20ζ − 101)2 b0 + (400ζ 2 − 8000ζ 3 + 8020ζ − 501)a40 + (−320000ζ 4 − 1172000ζ 3 + 400000ζ 2 + 1171950ζ − 94449)a30 + (−3200000ζ 5 − 7360000ζ 4 + 972000ζ 3 + 10459600ζ 2 + 2293530ζ − 3170800)a20 + (4800000ζ 5 − 36480000ζ 4 + 57612000ζ 3 + 39450000ζ 2 − 62415000ζ − 2954850)a0 − 1600000ζ 5 + 24160000ζ 4 + 60596000ζ 3 − 132200000ζ 2 − 58844500ζ + 108010000} . It took 0.015 second to compute the above shape basis in Maple 13 running on a 2.67 GHz PC with Intel Core i7 and 8 GB memory. It is observed that there are 3 polynomials and that the first polynomial, which is denoted by Sf (a0 ; ζ) hereafter, is 5th order in a0 and does not contain b1 and b0 . So, when the value of ζ is given, Sf (a0 ; ζ) is a univariate polynomial in a0 and computation of its real roots can be executed with guaranteed accuracy (Akritas, 1989). The second polynomial contains a0 , b1 and ζ and is in fact 1st order in b1 . The last polynomial contains a0 , b0 and ζ and is also 1st order in b0 . By equating each of the polynomials to 0, one can solve them for b1 and b0 , and the resulting expressions are polynomials in a0 and rational functions in ζ. When the value of ζ is given and a0 is computed, the values of b1 and b0 can be computed in a straightforward manner. In this way candidates of the optimal approximant are obtained. In fact singular cases occur when (ζ 2 − 1)(20ζ − 101)2 = 0, i.e., when ζ = ±1 or ζ = 101 20 ( = 5.05 ). The case ζ = −1 is excluded from the assumption that G(s) is stable. Also, when ζ = 101 20 , G(s) reduces to a 2nd order system: s + 150 G(s) = , (s + 1)(s + 100)
712
which should be treated separately (but is not considered further here). The only essential singularity thus happens when ζ = 1, where G(s) has multiple roots at s = −10: s + 150 G(s) = . (s + 1)(s + 10)2 This case is investigated first. The shape basis takes a slightly different form in this case (in addition to the obvious absence of ζ): {a40 + 179a30 + 70a20 − 150a0 − 1500 , 58320b1 − 151a30 − 26599a20 + 66420a0 + 77730 , 58320b0 − 2311a30 − 404599a20 + 1446660a0 − 1678350} . Notice that the first polynomial is now 4th order in a0 and that it has at most 4 real zeros. Of such zeros, positive ones are chosen to guarantee the stability of the approximants that are examined, and then the corresponding approximants are constructed. The H2 -norm of the error system is computed for each approximant and the one that yields the smallest H2 -norm error is found to be the optimal approximant. By computing real roots of the first polynomial, one gets −178.6036317 , 2.030690244 . From the only positive root, the optimal approximant is found to be 2.743100223 ˆ G(s) = , s + 2.030690244
ˆ = 0.3142081780. It is noted that achieving G(s) − G(s) 2 analysis of singular cases can systematically be dealt with by means of comprehensive Gr¨obner bases (Weispfenning, 1992), which is not pursued here. Now consider the generic case ζ 6= 1, 101 20 . When ζ is fixed, Sf has at most 5 real zeros and, in a similar manner, all the approximants calculated from the positive ones are to be examined. Let ζ = 51 = 0.2. Then, Sf (a0 ; ζ) reduces to 1 Sf a0 ; 5 = a50 + 157a40 + 708a30 − 3866a20 − 600a0 − 15000 , which has three real roots: −152.1808501 , −8.091702897 , 3.819877752 . The only positive root yields the best approximant 4.526338399 ˆ G(s) = s + 3.819877752
ˆ = 1.376827463. The Bode that achieves G(s) − G(s) 2 ˆ diagrams of G(s) and G(s) are shown in Fig. 1. Further examination of Sf (a0 ; ζ) gives more information on the optimal approximant. The constant term of Sf (a0 ; ζ) when seen as a polynomial in a0 is Sf (0; ζ) = −1500, which indicates that none of the roots of Sf (a0 ; ζ) crosses the origin as ζ changes. Also the discriminant of Sf (a0 ; ζ) has two roots 101 , ζ = −2.20557983679 , 20 and Sf (a0 ; 101 20 ) has a root a0 = −1 with multiplicity 2. Indeed the case ζ = 101 20 is excluded in the above analysis 101 and the cases 0 < ζ < 1, 1 < ζ < 101 20 and ζ > 20 are to be examined separately. However these facts still
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
20 G Ghat
0
Maginitude (dB)
Maginitude (dB)
20
−20 −40 −60 −2 10
−1
10
0
1
10 10 Frequency (rad/sec)
2
10
0 −20 −40 −60 −80 −2 10
3
10
−1
10
0
1
0
1
10 10 Frequency (rad/sec)
2
3
10
10
0
0
Phase (deg)
Phase (deg)
−90
−90
−180 −270 −360
−1
10
0
1
10 10 Frequency (rad/sec)
2
10
−450 −2 10
3
10
−1
10
10 10 Frequency (rad/sec)
2
3
10
10
ˆ cand,i (s) (ζ = 0.2) Fig. 2. Bode diagrams of G(s) and G
ˆ Fig. 1. Bode diagrams of G(s) and G(s) (ζ = 0.2) indicate that Sf (a0 ; ζ) has only one positive root when ζ > 0. Therefore, in order to find the optimal approximant of G(s), one has to compute the (only) positive root of Sf (a0 ; ζ) and substitute the value in the obtained expressions for other coefficients. It is noted that, as parametric expressions are available, the hardest part is the computation of the positive root of the 5th order polynomial Sf (a0 ; ζ) for a particular value of ζ.
0 Maginitude (dB)
−180 −2 10
G Ghat cand1 Ghat cand2, true optimal Ghat cand3
G Ghat
−20 −40 −60 −80 −2 10
−1
10
4.2 Another 3rd Order System
0
1
0
1
10 10 Frequency (rad/sec)
2
10
3
10
Seeing the previous example, one may wonder whether there will always be only one candidate for H2 model reduction from a 3rd order system to a 1st order one. The following example shows otherwise and in fact the number of candidates can vary as parameter values change: 50 1 − 2 , G(s) = s+1 s + 20ζs + 100 ˆ where ζ > 0. Write D(s) and Q(s) as in (2), and a set of polynomials similar to (3) is found. Again the Gr¨obner basis in shape form is computed for the generic case. Singular cases are when ζ = ±1, 101 20 , which are the same as in the previous example. The polynomial relating ζ and a0 is found to be Sf (a0 ; ζ) = a50 +(40ζ−151)a40 +(400ζ 2 −1040ζ−100)a30 + (−400ζ 2 + 2000ζ + 4650)a20 + (−5000ζ + 20000)a0 − 5000 in this case. This Sf (a0 ; ζ) is different from that in the previous example in that it can have more than one positive roots for some values of ζ. More specifically, when 0 < ζ < 2.769183502, Sf (a0 ; ζ) has 3 positive roots, meaning that there are 3 candidates. For instance, when ζ = 0.2, there are 3 positive roots 0.2472246004 , 6.571091706 , 144.7705454 , and 3 candidates derived from them are ˆ cand,1 (s) = 0.1517839124 , G s + 0.2472246004 ˆ cand,2 (s) = −2.141745375 , G s + 6.571091706 ˆ cand,3 (s) = 1.317210252 , G s + 144.7705454
713
Phase (deg)
0
−90
−180 −2 10
−1
10
10 10 Frequency (rad/sec)
2
10
3
10
ˆ Fig. 3. Bode diagrams of G(s) and G(s) (ζ = 4) respectively. The H2 -norm of the error system for each candidate is computed as
ˆ cand,1 (s) = 1.620501487 ,
G(s) − G
2 ˆ cand,2 (s) = 1.524330906 ,
G(s) − G
2
ˆ cand,3 (s) = 1.632980912 .
G(s) − G 2
These imply that the second candidate is the optimal one. ˆ cand,1 (s) and G ˆ cand,3 (s) try to The other two candidates G approximate the low frequency and high frequency ranges, respectively, while the true optimal one approximates both ranges reasonably well, still leaving out the peak around ω = 10 (Fig. 2). For 0.8635884203 < ζ < 2.769183502, ζ 6= 1, of 3 positive roots, the smallest one gives the optimal approximant, which mimics the low frequency range just as ˆ cand,1 (s) above. At ζ = 2.769183502, 2 of 3 positive G roots coincide and there are 2 distinct positive roots. For ζ > 2.769183502, there is only one positive root, which yields the global optimal approximant. Fig. 3 shows the case ζ = 4, where the optimal approximant derived from the only candidate is
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
40
−40 −60 −80 −2 10
Phase (deg)
Maginitude (dB)
−20
G Ghat1, true optimal Ghat2, true optimal −1
10
0
1
10 10 Frequency (rad/sec)
2
10
20 0 −20 −40 −60 −2 10
3
10
0
90
−90
0
Phase (deg)
Maginitude (dB)
0
−180 −270
−450 −2 10
−1
10
0
1
10 10 Frequency (rad/sec)
2
10
−270 −2 10
3
10
Fig. 4. Bode diagrams of G(s) and two optimal approxiˆ i (s) (ζ ≃ 0.8635884203) mants G 0.4011929392 ˆ . G(s) = s + 0.7929034808 Further investigation shows that the optimal approximant suddenly changes as ζ moves from 0 to 2.769183502. When ζ is close to 0, the second smallest positive root yields the optimal approximant. However candidates from the first and second smallest positive roots achieve the same error value when ζ takes the values of the smallest positive root of the following 10th order polynomial: 1366614474752ζ 10 + 21471716802560ζ 9 − 445896339292160ζ 8 + 2738751072378880ζ 7 − 7909799050741248ζ 6 + 10585295762298240ζ 5 − 3049368124144448ζ 4 − 5255439012035280ζ 3 − 125911367649880ζ 2 + 6934013226436280ζ − 3465965584594091 , which is approximately 0.8635884203. In other words there are two optimal approximants (Fig. 4): ˆ 1 (s) = 0.1730192650 , G s + 0.2878356873 ˆ 2 (s) = −0.8003214000 , G
s + 6.158636524
ˆ 1 (s) = G(s) − G ˆ 2 (s) = 0.5711503852 .
G(s) − G 2 2 This example shows that the H2 model reduction problem is a non-convex problem and further that some local optima can be a global one depending on the original system. 4.3 A 4th Order System Here the following 4th order system with a parameter ζ is approximated by a 2nd order system: 100 s2 + 4s + 1 · , G(s) = 2 s + 2ζs + 1 s2 + 2s + 100 where ζ > 0. The gain plot of G(s) has two peaks around ω = 1 and ω = 10 when 0 < ζ . 2 (Figs. 5 - 7). The height of the first peak around ω = 1 is determined by the value of ζ. It is anticipated that the 2nd order approximant would mimic one of the peaks and that the peak that is chosen
714
0
1
0
1
10 10 Frequency (rad/sec)
2
10
3
10
−90 −180
−360
−1
10
G Ghat cand1, true optimal Ghat cand2 Ghat cand3 −1
10
10 10 Frequency (rad/sec)
2
10
3
10
ˆ cand,i (s) (ζ = 0.05) Fig. 5. Bode diagrams of G(s) and G is determined by the height of the peak, or, equivalently, the value of ζ. In the same way as the 3rd order case in the previous subsections, a 17th order polynomial Sf (a0 ; ζ) in a0 is obtained, indicating that there are at most 17 candidates to be examined. The time required for computation of the shape basis was 136.064 seconds. 1 , the computation For comparison, when ζ is fixed to be 20 for this non-parametric case took 0.094 second. The optimal approximant is computed for the case of ζ = 1 20 = 0.05. The first peak is the higher in this case (Fig. 5). 1 ) has 8 positive roots. Then, It is found that Sf (a0 ; 20 a1 is computed, and there remains 3 stable candidates (candidates with a0 and a1 both being positive): 4.341427427s + 0.1013976740 ˆ cand,1 (s) = G , s2 + 0.1102306536s + 1.000448550 7.936404921s + 38.58538747 ˆ cand,2 (s) = G , s2 + 5.177755535s + 8.904076669 ˆ cand,3 (s) = −3.992932447s + 131.3864618 . G s2 + 2.835320555s + 97.36843591 The H2 -norm of the error system for each candidate is computed as
ˆ cand,1 (s) = 4.989215752 ,
G(s) − G
2 ˆ cand,2 (s) = 9.391652738 ,
G(s) − G
2
ˆ cand,3 (s) = 8.738207419 .
G(s) − G 2
ˆ cand,1 (s) is the optimal approximant This concludes that G and it can be seen that this approximant indeed shows a similar peak around ω = 1 in the gain plot (Fig. 5). Now the approximant for the case of ζ = 45 = 0.8 is sought. As the 17th order polynomial in a0 is already computed, one can reuse it with ζ = 45 . The main task is again the computation of positive roots of the polynomial Sf (a0 ; 54 ). There are 6 such roots in this case, and also there are 3 stable candidate approximants: 8.862740441s + 1.957116242 ˆ cand,1 (s) = G , s2 + 3.702075268s + 1.977288324 8.847950895s + 6.903893027 ˆ cand,2 (s) = G , s2 + 4.092974819s + 4.568928962 ˆ cand,3 (s) = −2.128943078s + 117.8865180 . G s2 + 2.354468593s + 98.67819435
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
5. CONCLUSION
Maginitude (dB)
40 20
This paper has considered the parametric H2 model reduction problem and solved it by means of Gr¨obner bases. The approach presented in this papers allows one to obtain the global optimal solution and also to analyse how the optimal approximant changes as parameter values change. While the computation times for the examples in this paper may be acceptable, the computational burden is still demanding. The development of more efficient computation approaches is thus anticipated.
0 −20 −40 −60 −2 10
−1
10
0
1
10 10 Frequency (rad/sec)
2
10
3
10
90
Phase (deg)
0
REFERENCES
−90 −180 −270 −2 10
G Ghat cand1 Ghat cand2 Ghat cand3, true optimal −1
10
0
1
10 10 Frequency (rad/sec)
2
10
3
10
ˆ cand,i (s) (ζ = 0.8) Fig. 6. Bode diagrams of G(s) and G
Maginitude (dB)
40 20 0 −20 −40 −60 −2 10
−1
10
0
1
0
1
10 10 Frequency (rad/sec)
2
10
3
10
Phase (deg)
90 0 −90 −180 −270 −2 10
G Ghat1, true optimal Ghat2, true optimal −1
10
10 10 Frequency (rad/sec)
2
10
3
10
Fig. 7. Bode diagrams of G(s) and two optimal approxiˆ i (s) (ζ ≃ 0.1389679763) mants G The H2 -norm of the error system for each candidate is computed as
ˆ cand,1 (s) = 4.627093455 ,
G(s) − G
2 ˆ cand,2 (s) = 4.630594568 ,
G(s) − G
2 ˆ cand,3 (s) = 1.187434431 .
G(s) − G 2
Since the second peak around ω = 10 is the higher in this ˆ cand,3 (s) that mimics this peak case, the approximant G achieves the optimum (Fig. 6). As may be predicted, for some value of ζ in between, there exist two global optimal approximants. In the case of ζ ≃ 0.1389679763, there are two approximants that achieve the global optimum (Fig. 7): 5.041464117s + 0.1512657375 ˆ 1 (s) = G , s2 + 0.3771111353s + 1.004477627 ˆ 2 (s) = −3.743921504s + 129.4998877 , G s2 + 2.755926095s + 97.56246240
ˆ 2 (s) = 4.932531593 . ˆ 1 (s) = G(s) − G
G(s) − G 2 2 Even in such a case, the approach can identify all the candidates and find the optimal approximant (up to the specified accuracy).
715
Akritas, A.G. (1989). Elements of Computer Algebra with Applications. John Wiley & Sons, New York, NY. Anderson, P.W. (1978). Local moments and localized states. Reviews of Modern Physics, 50(2), 191–201. Antoulas, A.C. (2005). Approximation of Large-scale Dynamical Systems. Advances in Design and Control. Society for Industrial and Applied Mathematics, Philadelphia. Cantoni, M. (2001). On model reduction in the ν-gap metric. In Proceedings of the 40th IEEE Conference on Decision and Control, 3665–3670. Orlando, Florida. Cox, D., Little, J., and O’Shea, D. (2007). Ideals, Varieties, and Algorithms. Springer, New York, NY, 3rd edition. Ghafoor, A. and Sreeram, V. (2008). A survey/review of frequency-weighted balanced model reduction techniques. Journal of Dynamic Systems, Measurement, and Control, 130(6). Glover, K. (1984). All optimal Hankel-norm approximations of linear multivariable systems and their L∞ -error bounds. International Journal of Control, 39(6), 1115– 1193. Gugercin, S., Antoulas, A.C., and Beattie, C. (2008). H2 model reduction for large-scale linear dynamical systems. SIAM Journal on Matrix Analysis and Applications, 30(2), 609–638. Hanzon, B. and Maciejowski, J.M. (1996). Constructive algebra methods for the L2 -problem for stable linear systems. Automatica, 32(12), 1645–1657. Kanno, M. (2010). H2 model reduction using an algebraic approach. In Proceedings of the SICE Annual Conference 2010, 2705–2708. Taipei, Taiwan. Meier III, L. and Luenberger, D.G. (1967). Approximation of linear constant systems. IEEE Transactions on Automatic Control, 12(5), 585–588. Moore, B.C. (1981). Principal component analysis in linear systems: Controllability, observability, and model reduction. IEEE Transactions on Automatic Control, 26(1), 17–32. Obinata, G. and Anderson, B.D.O. (2001). Model Reduction for Control System Design. Communications and Control Engineering. Springer-Verlag, London. Skelton, R.E., Iwasaki, T., and Grigoriadis, K. (1998). A Unified Algebraic Approach to Linear Control Design. Taylor & Francis Systems and Control Book Series. Taylor & Francis, London. Weispfenning, V. (1992). Comprehensive Gr¨obner bases. Journal of Symbolic Computation, 14(1), 1–29.
H2 Model Reduction by Means of Gr¨ obner Bases and its Extension to the Parametric Case Masaaki Kanno ∗ ∗
Institute of Science and Technology, Academic Assembly, Niigata University, Niigata, 950-2181, Japan (e-mail: [email protected]).
Abstract: In this paper a solution approach for parametric H2 model reduction based on Gr¨obner bases is presented. The approach utilizes a na¨ıve parametrization of the approximant, which leads to a rather simple set of algebraic equations that is solved by means of Gr¨obner bases. This allows one to deal with the parametric case, and analysis of the optimal solution in the presence of parameters becomes possible. Numerical examples show that the optimal approximant may change discontinuously as a parameter value varies, also showing the existence of systems with two optimal approximants. Keywords: H2 model reduction; algebraic methods; parametric case; multiple optimal solutions 1. INTRODUCTION Modelling is instrumental in modern science and engineering. A good model elucidates the nature of the object under investigation and helps one to construct a system that interacts with it. The meaning of a ‘good’ model can vary. An accurate model that reproduces the behaviour of the true object precisely may be beneficial in simulation, but can conceal beneath too much detail the real working mechanisms that want examination (Anderson, 1978). In modern control system design techniques, a detailed and complicated model may result in a controller of high orders, which can cause implementation issues due to computational complexity and numerical difficulties and may also be fragile against perturbation. The problem of model reduction has long been a subject of research in the systems and control field (Obinata and Anderson, 2001). In order to derive a simplified model from a complicated one, numerous model reduction techniques have been developed. One of popular methods is balanced truncation initiated by Moore (1981), which is further advanced by frequency-weighting (Ghafoor and Sreeram, 2008). Another popular approach is the Hankel norm approximation method (Glover, 1984). The advantage of those approaches is that they admit constructive algorithms. By setting some cost functions to measure the ‘closeness’ of the original system and the simplified one, model reduction can be formulated as optimization problems. Solution approaches appear for H2 and H∞ model reduction, e.g., Skelton et al. (1998). Another optimization problem setting relevant in feedback control is ν-gap model reduction investigated in Cantoni (2001) since the ν-gap is an appropriate metric to measure the closeness of two systems in a feedback loop. This paper considers H2 model reduction. The use of the H2 -norm to evaluate the fidelity of the reduced order model is advantageous in that this criterion admits both the time-domain and frequency-domain interpretations. 978-3-902823-06-9/12/$20.00 © 2012 IFAC
710
Investigation on this problem has a long history, and a wide range of techniques have been proposed. Many approaches rely on the first order necessary condition for optimality (Meier III and Luenberger, 1967) and its extension to multivariate case (Gugercin et al., 2008). A recent result numerically finds approximants for systems with thousands of states effectively (Antoulas, 2005; Gugercin et al., 2008). The problem of H2 model reduction is a nonconvex problem and it is intrinsically difficult to find the global optimal approximant. Even though a numerically effective approach can deal with a very large-scale system, the quality of the resulting approximant thus depends on the choice of the initial guess. The approach employed in this paper is based on an algebraic approach, more specifically Gr¨obner bases (Cox et al., 2007). The advantage of such approach over typical numerical approaches is that it may achieve global optimality and also can handle the parametric case. The drawback (and a major one) is computational complexity due to the reliance of symbolic manipulation and exact arithmetic. The aims of the paper is first to present the current state of an algebraic approach and further to investigate the H2 model reduction problem with the aid of such approach suggesting how difficult the problem can be. The paper is organized as follows. Section 2 states the problem formulation, followed by the suggested solution approach in Section 3. Three numerical examples are given in Section 4, where the approach is demonstrated and further it is shown how the number of candidates and the optimal reduced order model may change as a parameter value changes. Also systems with two global optimal solutions are presented. Some concluding remarks are made in Section 5. 2. PROBLEM FORMULATION The problem of H2 model reduction considered in this paper is formulated as follows. Given an n-th order stable 10.3182/20120711-3-BE-2027.00160
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
single-input, single-output system G(s), the task is to find ˆ an m-th order stable system G(s) (where m < n) that minimizes
ˆ ,
G(s) − G(s) 2 i.e., to find the best approximant in the H2 sense. The problem has been studied for decades and a wide variety of mathematical tools have been applied to solve the problem. An algebraic approach is utilized in this paper and one of early work in that direction is Hanzon and Maciejowski (1996). 3. SOLUTION APPROACH The solution approach which is employed here and initiated in Kanno (2010) relies on an algebraic method. The coefficients of the denominator of the approximant transfer function and those of an auxiliary polynomial are taken as unknown variables, and a set of algebraic equations in terms of those variables are derived from an optimality condition. Then a set of solutions, which constitutes a set of candidates, are calculated by means of an important algebraic method called Gr¨obner bases (Cox et al., 2007). Of such candidates, a stable system with real coefficients that achieves the minimal H2 -norm is the true optimal approximant. The main difference of the approach employed here from that of Hanzon and Maciejowski (1996) is a simple parametrization of the approximant. The approach by Hanzon and Maciejowski (1996) exploits the Schwarz-like canonical form which simplifies the stability criterion. The approach in this paper requires an extra effort and depends upon, e.g., the Routh-Hurwitz criterion to guarantee the stability of the approximant. Nevertheless simpler equations are obtained and the approach can deal with the parametric case with relative ease, which allows one to appreciate the strength of the algebraic approach. For instance, one may analyse how the optimal approximant changes as parameter values vary. Write the original system G(s) and its reduced order ˆ approximant G(s) to be sought as ˆ (s) N N (s) ˆ , G(s) = , G(s) = ˆ D(s) D(s) ˆ where D(s) and D(s) are stable monic polynomials of ˆ (s) are degree n and m, respectively, and N (s) and N polynomials of degree strictly less than n and m, respecˆ tively. For the optimal approximant G(s), there exists an auxiliary polynomial Q(s) of degree at most n−m−1 that satisfies the following equation (Meier III and Luenberger, 1967): 2 ˆ ˆ ˆ (s)D(s) . Q(s) = N (1) N (s)D(s) − D(−s) In this paper it is a standing assumption that the optimal ˆ reduced order model can be obtained by finding D(s) and Q(s) satisfying (1). Write ˆ D(s) = sm + am−1 sm−1 + · · · + a1 s + a0 , Q(s) = bn−m−1 sn−m−1 + · · · + b1 s + b0 , and compute the remainder R(s) on dividing the left hand side of (1) by D(s). Since R(s) has to be identically 0, its
711
coefficients, which are polynomials in ai ’s and bj ’s, have to be all zero. One thus obtains a set of algebraic equations. ˆ (s) can be computed Once the equations are solved, N by (1). In order to solve the obtained equations, a Gr¨obner basis of the polynomial parts of the equations is computed, which yields an equivalent set of equations that is easier to ‘solve’. It is further assumed that the system of equations has a finite number of solutions, i.e., the solution set is 0-dimensional. The Gr¨obner basis with respect to the lexicographic ordering bn−m−1 ≻ bn−m−2 ≻ · · · ≻ b0 ≻ am−1 ≻ am−2 ≻ · · · ≻ a0 is generically in the shape form: { Sf (a0 ) , a1 − h1 (a0 ) , . . . , am−1 − hm−1 (a0 ) , b0 − hm (a0 ) , . . . , bn−m−1 − hn−1 (a0 ) } . The set of zeros is then obtained by firstly finding the real roots of Sf (a0 ) and then substituting the values in the expressions hi (a0 ). Some of the zeros may yield unstable systems, which are to be discarded at this point. One then ˆ computes the H2 -norm of the error system G(s) − G(s) for each candidate and picks up the candidate giving the smallest H2 error, which is the optimal approximant. In the parametric case, Sf (a0 ) contains parameters, namely, Sf (a0 ) is a polynomial in a0 whose coefficients are functions in parameters. If parameters enter in polynomial form in the coefficients of G(s), then Sf (a0 ) will also be a polynomial in parameters. Once a parametric Sf (a0 ) is obtained, one can exploit such expression to compute the optimal approximant for each set of parameter values, and thus can utilize the result of off-line computation and avoid recomputation from the start. Namely one has to substitute the parameter with the value and compute positive roots of a polynomial. Then a finite number of candidates are obtained, and the H2 -norm of the error system for each case is to be calculated for determining the true optimal approximant. Further analysis may provide information on how many candidates can be obtained. This point is demonstrated in numerical examples in the subsequent section. 4. NUMERICAL EXAMPLES 4.1 A 3rd Order System The following 3rd order system with a parameter ζ is approximated by a 1st order system: 1 50 G(s) = + 2 s+1 s + 20ζs + 100 s2 + (20ζ + 50)s + 150 N (s) = 3 =: , 2 s + (20ζ + 1)s + (20ζ + 100)s + 100 D(s) where ζ > 0 (stability requirement). The denominaˆ ˆ tor D(s) of the approximant G(s) and the auxiliary polynomial Q(s) are written as ˆ D(s) = s + a0 , Q(s) = b1 s + b0 . (2) The quotient and the remainder on dividing the left hand side of (1) by D(s) are
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
1 − b1 , 2a0 b1 + a0 − b0 + (20ζ + 1)b1 + 49 s2 + −a20 b1 + 2a0 b0 + (20ζ + 50)a0 + (20ζ + 100)b1 −20ζ + 50 s − a20 b0 + 150a0 + 100b1 − 100 , respectively. The remainder has to be identically 0, and thus what is to be found is the zeros of the following ideal: {2a0 b1 + a0 − b0 + (20ζ + 1)b1 + 49 , − a20 b1 + 2a0 b0 + (20ζ + 50)a0 + (20ζ + 100)b1 − 20ζ + 50 , − a20 b0 + 150a0 + 100b1 − 100} . (3) Computing the Gr¨obner basis with respect to the lexicographic ordering b1 ≻ b0 ≻ a0 , one can get the Gr¨obner basis in shape form for the generic case: {a50 + (40ζ + 149)a40 + (400ζ 2 + 960ζ + 500)a30 + (−400ζ 2 + 6000ζ − 5050)a20 − 3000ζa0 − 15000 , 200(ζ 2 − 1)(20ζ − 101)2 b1 + (−400ζ 2 + 20ζ + 299)a40 + (−16000ζ 3 − 58600ζ 2 + 15950ζ + 44151)a30 + (−160000ζ 4 − 368000ζ 3 + 8600ζ 2 + 433530ζ + 80200)a20 + (240000ζ 4 − 1824000ζ 3 + 2910600ζ 2 − 3135450 + 1821000ζ)a0 − 80000ζ 4 + 1408000ζ 3 + 979800ζ 2 − 1256500ζ − 929800 , 200(ζ 2 − 1)(20ζ − 101)2 b0 + (400ζ 2 − 8000ζ 3 + 8020ζ − 501)a40 + (−320000ζ 4 − 1172000ζ 3 + 400000ζ 2 + 1171950ζ − 94449)a30 + (−3200000ζ 5 − 7360000ζ 4 + 972000ζ 3 + 10459600ζ 2 + 2293530ζ − 3170800)a20 + (4800000ζ 5 − 36480000ζ 4 + 57612000ζ 3 + 39450000ζ 2 − 62415000ζ − 2954850)a0 − 1600000ζ 5 + 24160000ζ 4 + 60596000ζ 3 − 132200000ζ 2 − 58844500ζ + 108010000} . It took 0.015 second to compute the above shape basis in Maple 13 running on a 2.67 GHz PC with Intel Core i7 and 8 GB memory. It is observed that there are 3 polynomials and that the first polynomial, which is denoted by Sf (a0 ; ζ) hereafter, is 5th order in a0 and does not contain b1 and b0 . So, when the value of ζ is given, Sf (a0 ; ζ) is a univariate polynomial in a0 and computation of its real roots can be executed with guaranteed accuracy (Akritas, 1989). The second polynomial contains a0 , b1 and ζ and is in fact 1st order in b1 . The last polynomial contains a0 , b0 and ζ and is also 1st order in b0 . By equating each of the polynomials to 0, one can solve them for b1 and b0 , and the resulting expressions are polynomials in a0 and rational functions in ζ. When the value of ζ is given and a0 is computed, the values of b1 and b0 can be computed in a straightforward manner. In this way candidates of the optimal approximant are obtained. In fact singular cases occur when (ζ 2 − 1)(20ζ − 101)2 = 0, i.e., when ζ = ±1 or ζ = 101 20 ( = 5.05 ). The case ζ = −1 is excluded from the assumption that G(s) is stable. Also, when ζ = 101 20 , G(s) reduces to a 2nd order system: s + 150 G(s) = , (s + 1)(s + 100)
712
which should be treated separately (but is not considered further here). The only essential singularity thus happens when ζ = 1, where G(s) has multiple roots at s = −10: s + 150 G(s) = . (s + 1)(s + 10)2 This case is investigated first. The shape basis takes a slightly different form in this case (in addition to the obvious absence of ζ): {a40 + 179a30 + 70a20 − 150a0 − 1500 , 58320b1 − 151a30 − 26599a20 + 66420a0 + 77730 , 58320b0 − 2311a30 − 404599a20 + 1446660a0 − 1678350} . Notice that the first polynomial is now 4th order in a0 and that it has at most 4 real zeros. Of such zeros, positive ones are chosen to guarantee the stability of the approximants that are examined, and then the corresponding approximants are constructed. The H2 -norm of the error system is computed for each approximant and the one that yields the smallest H2 -norm error is found to be the optimal approximant. By computing real roots of the first polynomial, one gets −178.6036317 , 2.030690244 . From the only positive root, the optimal approximant is found to be 2.743100223 ˆ G(s) = , s + 2.030690244
ˆ = 0.3142081780. It is noted that achieving G(s) − G(s) 2 analysis of singular cases can systematically be dealt with by means of comprehensive Gr¨obner bases (Weispfenning, 1992), which is not pursued here. Now consider the generic case ζ 6= 1, 101 20 . When ζ is fixed, Sf has at most 5 real zeros and, in a similar manner, all the approximants calculated from the positive ones are to be examined. Let ζ = 51 = 0.2. Then, Sf (a0 ; ζ) reduces to 1 Sf a0 ; 5 = a50 + 157a40 + 708a30 − 3866a20 − 600a0 − 15000 , which has three real roots: −152.1808501 , −8.091702897 , 3.819877752 . The only positive root yields the best approximant 4.526338399 ˆ G(s) = s + 3.819877752
ˆ = 1.376827463. The Bode that achieves G(s) − G(s) 2 ˆ diagrams of G(s) and G(s) are shown in Fig. 1. Further examination of Sf (a0 ; ζ) gives more information on the optimal approximant. The constant term of Sf (a0 ; ζ) when seen as a polynomial in a0 is Sf (0; ζ) = −1500, which indicates that none of the roots of Sf (a0 ; ζ) crosses the origin as ζ changes. Also the discriminant of Sf (a0 ; ζ) has two roots 101 , ζ = −2.20557983679 , 20 and Sf (a0 ; 101 20 ) has a root a0 = −1 with multiplicity 2. Indeed the case ζ = 101 20 is excluded in the above analysis 101 and the cases 0 < ζ < 1, 1 < ζ < 101 20 and ζ > 20 are to be examined separately. However these facts still
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
20 G Ghat
0
Maginitude (dB)
Maginitude (dB)
20
−20 −40 −60 −2 10
−1
10
0
1
10 10 Frequency (rad/sec)
2
10
0 −20 −40 −60 −80 −2 10
3
10
−1
10
0
1
0
1
10 10 Frequency (rad/sec)
2
3
10
10
0
0
Phase (deg)
Phase (deg)
−90
−90
−180 −270 −360
−1
10
0
1
10 10 Frequency (rad/sec)
2
10
−450 −2 10
3
10
−1
10
10 10 Frequency (rad/sec)
2
3
10
10
ˆ cand,i (s) (ζ = 0.2) Fig. 2. Bode diagrams of G(s) and G
ˆ Fig. 1. Bode diagrams of G(s) and G(s) (ζ = 0.2) indicate that Sf (a0 ; ζ) has only one positive root when ζ > 0. Therefore, in order to find the optimal approximant of G(s), one has to compute the (only) positive root of Sf (a0 ; ζ) and substitute the value in the obtained expressions for other coefficients. It is noted that, as parametric expressions are available, the hardest part is the computation of the positive root of the 5th order polynomial Sf (a0 ; ζ) for a particular value of ζ.
0 Maginitude (dB)
−180 −2 10
G Ghat cand1 Ghat cand2, true optimal Ghat cand3
G Ghat
−20 −40 −60 −80 −2 10
−1
10
4.2 Another 3rd Order System
0
1
0
1
10 10 Frequency (rad/sec)
2
10
3
10
Seeing the previous example, one may wonder whether there will always be only one candidate for H2 model reduction from a 3rd order system to a 1st order one. The following example shows otherwise and in fact the number of candidates can vary as parameter values change: 50 1 − 2 , G(s) = s+1 s + 20ζs + 100 ˆ where ζ > 0. Write D(s) and Q(s) as in (2), and a set of polynomials similar to (3) is found. Again the Gr¨obner basis in shape form is computed for the generic case. Singular cases are when ζ = ±1, 101 20 , which are the same as in the previous example. The polynomial relating ζ and a0 is found to be Sf (a0 ; ζ) = a50 +(40ζ−151)a40 +(400ζ 2 −1040ζ−100)a30 + (−400ζ 2 + 2000ζ + 4650)a20 + (−5000ζ + 20000)a0 − 5000 in this case. This Sf (a0 ; ζ) is different from that in the previous example in that it can have more than one positive roots for some values of ζ. More specifically, when 0 < ζ < 2.769183502, Sf (a0 ; ζ) has 3 positive roots, meaning that there are 3 candidates. For instance, when ζ = 0.2, there are 3 positive roots 0.2472246004 , 6.571091706 , 144.7705454 , and 3 candidates derived from them are ˆ cand,1 (s) = 0.1517839124 , G s + 0.2472246004 ˆ cand,2 (s) = −2.141745375 , G s + 6.571091706 ˆ cand,3 (s) = 1.317210252 , G s + 144.7705454
713
Phase (deg)
0
−90
−180 −2 10
−1
10
10 10 Frequency (rad/sec)
2
10
3
10
ˆ Fig. 3. Bode diagrams of G(s) and G(s) (ζ = 4) respectively. The H2 -norm of the error system for each candidate is computed as
ˆ cand,1 (s) = 1.620501487 ,
G(s) − G
2 ˆ cand,2 (s) = 1.524330906 ,
G(s) − G
2
ˆ cand,3 (s) = 1.632980912 .
G(s) − G 2
These imply that the second candidate is the optimal one. ˆ cand,1 (s) and G ˆ cand,3 (s) try to The other two candidates G approximate the low frequency and high frequency ranges, respectively, while the true optimal one approximates both ranges reasonably well, still leaving out the peak around ω = 10 (Fig. 2). For 0.8635884203 < ζ < 2.769183502, ζ 6= 1, of 3 positive roots, the smallest one gives the optimal approximant, which mimics the low frequency range just as ˆ cand,1 (s) above. At ζ = 2.769183502, 2 of 3 positive G roots coincide and there are 2 distinct positive roots. For ζ > 2.769183502, there is only one positive root, which yields the global optimal approximant. Fig. 3 shows the case ζ = 4, where the optimal approximant derived from the only candidate is
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
40
−40 −60 −80 −2 10
Phase (deg)
Maginitude (dB)
−20
G Ghat1, true optimal Ghat2, true optimal −1
10
0
1
10 10 Frequency (rad/sec)
2
10
20 0 −20 −40 −60 −2 10
3
10
0
90
−90
0
Phase (deg)
Maginitude (dB)
0
−180 −270
−450 −2 10
−1
10
0
1
10 10 Frequency (rad/sec)
2
10
−270 −2 10
3
10
Fig. 4. Bode diagrams of G(s) and two optimal approxiˆ i (s) (ζ ≃ 0.8635884203) mants G 0.4011929392 ˆ . G(s) = s + 0.7929034808 Further investigation shows that the optimal approximant suddenly changes as ζ moves from 0 to 2.769183502. When ζ is close to 0, the second smallest positive root yields the optimal approximant. However candidates from the first and second smallest positive roots achieve the same error value when ζ takes the values of the smallest positive root of the following 10th order polynomial: 1366614474752ζ 10 + 21471716802560ζ 9 − 445896339292160ζ 8 + 2738751072378880ζ 7 − 7909799050741248ζ 6 + 10585295762298240ζ 5 − 3049368124144448ζ 4 − 5255439012035280ζ 3 − 125911367649880ζ 2 + 6934013226436280ζ − 3465965584594091 , which is approximately 0.8635884203. In other words there are two optimal approximants (Fig. 4): ˆ 1 (s) = 0.1730192650 , G s + 0.2878356873 ˆ 2 (s) = −0.8003214000 , G
s + 6.158636524
ˆ 1 (s) = G(s) − G ˆ 2 (s) = 0.5711503852 .
G(s) − G 2 2 This example shows that the H2 model reduction problem is a non-convex problem and further that some local optima can be a global one depending on the original system. 4.3 A 4th Order System Here the following 4th order system with a parameter ζ is approximated by a 2nd order system: 100 s2 + 4s + 1 · , G(s) = 2 s + 2ζs + 1 s2 + 2s + 100 where ζ > 0. The gain plot of G(s) has two peaks around ω = 1 and ω = 10 when 0 < ζ . 2 (Figs. 5 - 7). The height of the first peak around ω = 1 is determined by the value of ζ. It is anticipated that the 2nd order approximant would mimic one of the peaks and that the peak that is chosen
714
0
1
0
1
10 10 Frequency (rad/sec)
2
10
3
10
−90 −180
−360
−1
10
G Ghat cand1, true optimal Ghat cand2 Ghat cand3 −1
10
10 10 Frequency (rad/sec)
2
10
3
10
ˆ cand,i (s) (ζ = 0.05) Fig. 5. Bode diagrams of G(s) and G is determined by the height of the peak, or, equivalently, the value of ζ. In the same way as the 3rd order case in the previous subsections, a 17th order polynomial Sf (a0 ; ζ) in a0 is obtained, indicating that there are at most 17 candidates to be examined. The time required for computation of the shape basis was 136.064 seconds. 1 , the computation For comparison, when ζ is fixed to be 20 for this non-parametric case took 0.094 second. The optimal approximant is computed for the case of ζ = 1 20 = 0.05. The first peak is the higher in this case (Fig. 5). 1 ) has 8 positive roots. Then, It is found that Sf (a0 ; 20 a1 is computed, and there remains 3 stable candidates (candidates with a0 and a1 both being positive): 4.341427427s + 0.1013976740 ˆ cand,1 (s) = G , s2 + 0.1102306536s + 1.000448550 7.936404921s + 38.58538747 ˆ cand,2 (s) = G , s2 + 5.177755535s + 8.904076669 ˆ cand,3 (s) = −3.992932447s + 131.3864618 . G s2 + 2.835320555s + 97.36843591 The H2 -norm of the error system for each candidate is computed as
ˆ cand,1 (s) = 4.989215752 ,
G(s) − G
2 ˆ cand,2 (s) = 9.391652738 ,
G(s) − G
2
ˆ cand,3 (s) = 8.738207419 .
G(s) − G 2
ˆ cand,1 (s) is the optimal approximant This concludes that G and it can be seen that this approximant indeed shows a similar peak around ω = 1 in the gain plot (Fig. 5). Now the approximant for the case of ζ = 45 = 0.8 is sought. As the 17th order polynomial in a0 is already computed, one can reuse it with ζ = 45 . The main task is again the computation of positive roots of the polynomial Sf (a0 ; 54 ). There are 6 such roots in this case, and also there are 3 stable candidate approximants: 8.862740441s + 1.957116242 ˆ cand,1 (s) = G , s2 + 3.702075268s + 1.977288324 8.847950895s + 6.903893027 ˆ cand,2 (s) = G , s2 + 4.092974819s + 4.568928962 ˆ cand,3 (s) = −2.128943078s + 117.8865180 . G s2 + 2.354468593s + 98.67819435
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
5. CONCLUSION
Maginitude (dB)
40 20
This paper has considered the parametric H2 model reduction problem and solved it by means of Gr¨obner bases. The approach presented in this papers allows one to obtain the global optimal solution and also to analyse how the optimal approximant changes as parameter values change. While the computation times for the examples in this paper may be acceptable, the computational burden is still demanding. The development of more efficient computation approaches is thus anticipated.
0 −20 −40 −60 −2 10
−1
10
0
1
10 10 Frequency (rad/sec)
2
10
3
10
90
Phase (deg)
0
REFERENCES
−90 −180 −270 −2 10
G Ghat cand1 Ghat cand2 Ghat cand3, true optimal −1
10
0
1
10 10 Frequency (rad/sec)
2
10
3
10
ˆ cand,i (s) (ζ = 0.8) Fig. 6. Bode diagrams of G(s) and G
Maginitude (dB)
40 20 0 −20 −40 −60 −2 10
−1
10
0
1
0
1
10 10 Frequency (rad/sec)
2
10
3
10
Phase (deg)
90 0 −90 −180 −270 −2 10
G Ghat1, true optimal Ghat2, true optimal −1
10
10 10 Frequency (rad/sec)
2
10
3
10
Fig. 7. Bode diagrams of G(s) and two optimal approxiˆ i (s) (ζ ≃ 0.1389679763) mants G The H2 -norm of the error system for each candidate is computed as
ˆ cand,1 (s) = 4.627093455 ,
G(s) − G
2 ˆ cand,2 (s) = 4.630594568 ,
G(s) − G
2 ˆ cand,3 (s) = 1.187434431 .
G(s) − G 2
Since the second peak around ω = 10 is the higher in this ˆ cand,3 (s) that mimics this peak case, the approximant G achieves the optimum (Fig. 6). As may be predicted, for some value of ζ in between, there exist two global optimal approximants. In the case of ζ ≃ 0.1389679763, there are two approximants that achieve the global optimum (Fig. 7): 5.041464117s + 0.1512657375 ˆ 1 (s) = G , s2 + 0.3771111353s + 1.004477627 ˆ 2 (s) = −3.743921504s + 129.4998877 , G s2 + 2.755926095s + 97.56246240
ˆ 2 (s) = 4.932531593 . ˆ 1 (s) = G(s) − G
G(s) − G 2 2 Even in such a case, the approach can identify all the candidates and find the optimal approximant (up to the specified accuracy).
715
Akritas, A.G. (1989). Elements of Computer Algebra with Applications. John Wiley & Sons, New York, NY. Anderson, P.W. (1978). Local moments and localized states. Reviews of Modern Physics, 50(2), 191–201. Antoulas, A.C. (2005). Approximation of Large-scale Dynamical Systems. Advances in Design and Control. Society for Industrial and Applied Mathematics, Philadelphia. Cantoni, M. (2001). On model reduction in the ν-gap metric. In Proceedings of the 40th IEEE Conference on Decision and Control, 3665–3670. Orlando, Florida. Cox, D., Little, J., and O’Shea, D. (2007). Ideals, Varieties, and Algorithms. Springer, New York, NY, 3rd edition. Ghafoor, A. and Sreeram, V. (2008). A survey/review of frequency-weighted balanced model reduction techniques. Journal of Dynamic Systems, Measurement, and Control, 130(6). Glover, K. (1984). All optimal Hankel-norm approximations of linear multivariable systems and their L∞ -error bounds. International Journal of Control, 39(6), 1115– 1193. Gugercin, S., Antoulas, A.C., and Beattie, C. (2008). H2 model reduction for large-scale linear dynamical systems. SIAM Journal on Matrix Analysis and Applications, 30(2), 609–638. Hanzon, B. and Maciejowski, J.M. (1996). Constructive algebra methods for the L2 -problem for stable linear systems. Automatica, 32(12), 1645–1657. Kanno, M. (2010). H2 model reduction using an algebraic approach. In Proceedings of the SICE Annual Conference 2010, 2705–2708. Taipei, Taiwan. Meier III, L. and Luenberger, D.G. (1967). Approximation of linear constant systems. IEEE Transactions on Automatic Control, 12(5), 585–588. Moore, B.C. (1981). Principal component analysis in linear systems: Controllability, observability, and model reduction. IEEE Transactions on Automatic Control, 26(1), 17–32. Obinata, G. and Anderson, B.D.O. (2001). Model Reduction for Control System Design. Communications and Control Engineering. Springer-Verlag, London. Skelton, R.E., Iwasaki, T., and Grigoriadis, K. (1998). A Unified Algebraic Approach to Linear Control Design. Taylor & Francis Systems and Control Book Series. Taylor & Francis, London. Weispfenning, V. (1992). Comprehensive Gr¨obner bases. Journal of Symbolic Computation, 14(1), 1–29.