- Email: [email protected]
Controller reduction via minimum rank matrix approximation
Automatica 48 (2012) 1069–1076
Contents lists available at SciVerse ScienceDirect
Automatica journal homepage: www.elsevier.com/locate/automatica
Controller reduction via minimum rank matrix approximation✩ Kin Cheong Sou a,1 , Anders Rantzer b a
Automatic Control Lab and ACCESS Linnaeus Center, School of Electrical Engineering, Royal Institute of Technology, Osquldas väg 10, floor 6, SE-100 44 Stockholm, Sweden
b
Department of Automatic Control, Lund University, Box 118, SE-221 00 Lund, Sweden
article
abstract
info
Article history: Received 7 April 2011 Received in revised form 17 August 2011 Accepted 27 October 2011 Available online 30 March 2012
In this paper a controller reduction method for discrete-time linear time-invariant systems is described. Using the bounded-real lemma, the proposed method generates reduced controllers with closed loop stability and H∞ norm performance guarantee. Information of the full controller is used as a basis for reduction using singular value decomposition. This is different from traditional model reduction schemes such as weighted balanced truncation. Numerical assessment of the proposed method is given in the end of the paper. © 2012 Elsevier Ltd. All rights reserved.
Keywords: Model reduction Singular value decomposition
1. Introduction In this paper a controller reduction problem for discrete-time linear time-invariant systems depicted in Fig. 1 is considered. In this reduction problem, the feedback structure should remain and only the controller can be simplified. This can be viewed as a special case of the more general problem of structured model reduction (e.g. Sandberg (2010)). However, the controller reduction problem, by itself, has also been extensively studied because of its potential values. The main idea of the proposed controller reduction method can be explained as follows. Standard H∞ and H2 synthesis procedures involve solution of the following linear matrix inequality (LMI):
X 0
0 γI
A ≻ cl Ccl
Bcl Dcl
T
X
0
0 1 Acl Ccl I
γ
Bcl Dcl
(1)
where Acl , Bcl , Ccl and Dcl are affine functions of an unknown matrix L describing the controller, and X is another unknown as the ‘‘KYP certificate’’. The optimal L and X found by solving (1) often lead to high order controllers. We are interested in finding good approximations of low order. The approach is to use the same X as
✩ The material in this paper was partially presented at the 2010 American Control Conference, June 30–July 2, 2010, Baltimore, Maryland, USA and the 49th IEEE Conference on Decision and Control (CDC), December 15–17, 2010, Atlanta, Georgia, USA. This paper was recommended for publication in revised form by Associate Editor Richard D. Braatz under the direction of Editor Frank Allgöwer. E-mail addresses: [email protected] (K.C. Sou), [email protected] (A. Rantzer). 1 Tel.: +46 8 7907427; fax: +46 8 7907329.
0005-1098/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.automatica.2012.03.011
for the optimal high order controller, then replace L by a matrix Lˆ of minimal rank, subject to the inequality (1). Rank minimization is not equivalent to controller order reduction, but low rank leads to low order controllers and works well in examples. The above idea can be viewed as the middle ground between two existing groups of controller reduction methods, balancing computation effort and approximation quality. One group is balanced truncation/Hankel norm reduction and their weighted variants (e.g. Enns (1984), Anderson and Liu (1989), Goddard and Glover (1993), Zhou (1993), and Obinata and Anderson (2001)). These methods utilize full controller information. They are efficient and produce good local approximation results. The other group of known reduction methods is optimization based, typically formulated using the KYP lemma (e.g. Gahinet and Apkarian (1994), Iwasaki and Skelton (1994), Ankelhed, Helmersson, and Hansson (2006), Ajala, Schuler, and Allgower (2010), and Iwasaki and Skelton (1995a,b)). These methods involve more expensive computations, but a globally optimal reduced controller can be found if the corresponding optimization problems are solved to optimality. Full controller information is not always necessary in these optimization based methods. The purpose of this paper is not to present yet another attempt to solve the KYP problem to optimality. Rather, it describes a perspective on how to utilize the full controller information to obtain a reduced controller. This perspective is different from weighted balanced truncation, and it attempts to combine the benefits of weighted balanced truncation and KYP optimization methods. In addition, via numerical examples this paper provides some assessment of how well the investigated perspective can work in practice. Nevertheless, the full theoretical characterization of the limit of performance remains unknown in this paper.
1070
K.C. Sou, A. Rantzer / Automatica 48 (2012) 1069–1076 w
w
z P
u
y
z P
u
K
With the matrices defined in (4), the state space matrices of the reduced closed loop system P ⋆ Kˆ can be described as affine functions of the reduced controller Lˆ as
y
ˆ cl , A¯ + B¯ 2 Lˆ C¯ 2 A
K
Fig. 1. The block diagram on the left is the feedback interconnection of a plant P and a controller K considered in this paper. The problem of controller reduction is to find a lower order controller Kˆ to replace K , so that the performances of the reduced closed loop system (on the right) do not degrade too much.
The rest of the paper is organized as follows. In Section 2 the background material and notations are introduced. In Section 3 the proposed controller reduction problem is motivated in detail and then formulated. Then Section 4 describes the SVD based solution to the proposed problem. The order of the reduced controller turns out to be dependent on a ‘‘KYP matrix’’, the search of which is to be discussed in Section 5. Numerical assessments of the proposed method are presented in Section 6. 2. Background and notation 2.1. Controller reduction problem In this paper, all the systems involved (except some in Section 6) are discrete-time linear time-invariant systems. The main object of interest is the feedback setup of a plant P and a controller K in Fig. 1. The feedback closed loop system is denoted as P ⋆ K , and it is assumed to be stable. In the controller reduction problem in this paper, P and K are given data. A reduced controller Kˆ is sought as the minimum order controller whose corresponding closed loop system, denoted as P ⋆ Kˆ , has H∞ norm less than a userspecified upper bound. The controller reduction problem as stated is computationally difficult to solve, and therefore an approximate solution will be found by solving a different but related problem proposed in Section 3. 2.2. State space characterization of closed loop system In Fig. 1 the dimensions of the signals w, u, z , y are nw , nu , nz , ny respectively. The plant P is described in state space matrices as follows.
C1 P = (zI − A)−1 B1 C2
B2
D + 11 D21
D12 D22
(2)
where A ∈ Rn×n , D12 ∈ Rnz ×nu , D21 ∈ Rny ×nw and z ∈ C. In this paper, a technical assumption is made that D22 = 0 (see Doyle, Glover, Khargonekar, and Francis (1989) for detail). The given full controller K is described by its state space matrices as K = Ck (zI − Ak )−1 Bk + Dk . Similarly, the reduced ˆ k . In above, Ak ∈ Rnk ×nk , Dk ∈ controller Kˆ is Cˆ k (zI − Aˆ k )−1 Bˆ k + D ˆ nu ×ny n ׈ n n × ˆ k ∈ R u ny . The subsequent exposition R and Aˆ k ∈ R k k , D requires the definitions of the following two matrices.
Ak L, Ck
Bk Dk
Aˆ k ˆ and L , Cˆ k
Bˆ k ˆk . D
(3)
Subsequently, L and Lˆ are referred to as the full and reduced controller matrices, respectively. Also, define the data matrices (A¯ , B¯ 1 , B¯ 2 , etc.) A 0 B1 0 B2 0 0 0 Inˆ k 0 B¯ 1 B¯ 2 A¯ nˆ k ¯ ¯ 12 , C1 C1 D11 D 0 D11 0 D12 . ¯ 21 ¯ 22 C¯ 2 D D 0 Inˆ 0 0 0 C2
k
0
D21
0
0
(4)
¯ 12 Lˆ C¯ 2 Cˆ cl , C¯ 1 + D
ˆ cl , B¯ 1 + B¯ 2 Lˆ D¯ 21 B
(5)
ˆ cl , D11 + D¯ 12 Lˆ D¯ 21 . D
The dimensions of the closed loop state space matrices above depend on the dimension of Lˆ (i.e. nˆ k ). Also, let Acl , Bcl , Ccl and Dcl denote the state space matrices of the full closed loop system P ⋆ K . In the case where nˆ k = nk , which is the main focus of this paper, these matrices can also be written as in (5), with L replacing ˆ Finally, ∥P ⋆ K ∥∞ and ∥P ⋆ Kˆ ∥∞ denote the H∞ norms of the L. respective closed loop systems. 3. Basic problem formulation 3.1. Closed loop performance objectives Two performance objectives should be guaranteed by the reduced controller Kˆ . First, P ⋆ Kˆ should be stable. Second, for any given γ > ∥P ⋆ K ∥∞ , it holds that ∥P ⋆ Kˆ ∥∞ < γ (the lower bound on γ ensures that K can be one of the candidate reduced controllers). These two performance objectives can be specified using the state space matrices of P ⋆ Kˆ (which, according to (5), ˆ via the bounded-real lemma as follows. are functions of L), Theorem 1 (Zhou, Doyle, & Glover, 1995). Let γ be a positive scalar. In addition, denote the closed loop system corresponding to Lˆ (cf. (3)
ˆ cl and (5)) as P ⋆ Kˆ = Cˆ cl zI − A
−1
ˆ cl + Dˆ cl . The following two B
statements are equivalent.
• P ⋆ Kˆ is stable and ∥P ⋆ Kˆ ∥∞ < γ • there exists X ≻ 0 such that 0 ˆ ˆ cl Bˆ cl T X X 0 A 1 Acl ≻ ˆ ˆ cl 0 γI I 0 Ccl D Cˆ cl γ
(6a)
ˆ cl B ˆ cl . D
(6b)
Remark 1. The results in this paper rely on the structure in (6b). Since the continuous-time bounded-real lemma does not have a statement similar to (6b) in structure, this paper considers the discrete-time case only. 3.2. A convex restriction of controller matrix candidates It is a well-known difficulty to apply (6b) to simultaneously solve for the bound-real matrices X and reduced controller ˆ Refs. (Gahinet & Apkarian, 1994; Iwasaki & Skelton, matrices L. 1994) report a convex equivalence of (6b), but it remains nonconvex if an order constraint is imposed (i.e., nˆ k < nk ). Instead, this paper proposes to use the full controller matrix L to separate ˆ First X is found as a solution to the following the search of X and L. LMI
X 0
0 γI
A ≻ cl Ccl
Bcl Dcl
T
X
0
0 1 Acl Ccl I
γ
Bcl Dcl
(7)
where Acl , Bcl , Ccl and Dcl are defined in the discussion following (5). (7) ensures that X is admissible in the sense that at least L satisfies (6b) with X defined above. Once X is determined, the corresponding Lˆ is found by solving (6b). Imposing (7) leads to a restriction on the choice of X as opposed to letting it be an unknown in (6b). Nevertheless, often the full controller is
K.C. Sou, A. Rantzer / Automatica 48 (2012) 1069–1076
unnecessarily high order. Hence, it is sensible to expect some lower order controllers would also satisfy (6b) with X from (7) with an appropriate γ which leads to acceptable closed loop performances. In addition, imposing (7) is in line with local methods such as weighted balanced truncation — to utilize the full controller information as well as possible. In this paper, a detailed evaluation through numerical examples is conducted to assess the limit of performance imposed by (7), though a full theoretical investigation remains difficult. 3.3. Order specification via rank of the controller matrix Apart from the restrictiveness, (7) also complicates the search of controllers of reduced orders. In particular, the combination of (5) and (7) implies that nˆ k = nk (i.e. the A matrices of the full and reduced controllers have the same dimension). Hence, to achieve order reduction it would be desirable to directly ˆ However, this renders minimize order(Kˆ ) (which is a function of L). ˆ the search of L intractable. Instead, this paper proposes that rank(Lˆ ) should be minimized. This is based on the following result preliminarily from Sou and Rantzer (2010a,b) (now proved in the Appendix). Proposition 2. Let Kˆ be the reduced controller associated with a reduced controller matrix Lˆ as in (3), then
order(Kˆ ) ≤ min rank Aˆ k Aˆ k Cˆ k
≤ rank
Bˆ k ˆk D
Bˆ k
Aˆ k , rank Cˆ k
The reduced controller is found via the following steps: (1) Solve the LMI’s in (8c) to obtain X . The solutions to (8c) are typically non-unique, and there are different approaches to find the different X . The details are given in Section 5. (2) Once X is determined from (8c), solve the constrained rank minimization problem of (8a) and (8b). This problem can be solved efficiently using SVD. See Section 4. (3) Once Lˆ of the minimum rank is determined, form a minimal state space model from Lˆ as the proposed reduced order controller in the problem in Fig. 1. 4. Minimum rank controller matrix Lˆ The main issue in this section is to solve the constrained minimum rank problem (8a) and (8b), assuming that the matrix X has been determined from (8c). The main result here is that the seemingly intractable rank minimization problem in (8a) and (8b) can be shown to be equivalent to a SVD solvable classical matrix approximation problem to be defined in (16). To show the equivalence, first the matrix inequality in (8b) is simplified using Theorem 1 in Sou and Rantzer (2010c). Then, (8a) and the simplified version of (8b) becomes a classical minimum rank matrix approximation problem, solvable using SVD. 4.1. Simplification of matrix inequality The main objective of this subsection is to show that (8b) and a simplified constraint to be defined in (13) are equivalent. The development is as follows. Using on the (5), the last matrix
= rank(Lˆ ).
ˆ for the low Nevertheless, (7) might exclude the forming of such L’s order controllers. This remains an open question, but Section 6 provides an experimental indication that rank minimization need not be restrictive. 3.4. Formulation of the basic problem The data of the proposed problem are the plant matrices in (4) and L in (3) describing the full order controller. There is a user-specified scalar parameter γ such that γ > ∥P ⋆ K ∥∞ and the proposed problem would enforce that γ > ∥P ⋆ Kˆ ∥∞ , thus controlling the tradeoff between performance loss and controller order reduction. The reduced controller matrix Lˆ is the decision variable. Lˆ and L have the same dimensions (i.e. nˆ k = nk ), as discussed before. The proposed problem, in its basic form, can be formulated as minimize rank(Lˆ ) subject to
(8a)
Lˆ
X 0
0
γI
ˆ A ≻ ˆ cl Ccl
ˆ cl B ˆ cl D
T
X
0
0 ˆ 1 Acl I Cˆ cl
γ
ˆ cl B ˆ cl D
(8b)
where X ≻ 0 can be computed from the data by solving
X 0
0 γI
A ≻ cl Ccl
Bcl Dcl
T
X
0
0 1 Acl Ccl I
γ
1
(8c)
The definitions of the closed loop system matrices in (8b) and (8c) can be found in (5). Also, see Sections 3.1–3.3 for the explanations of (8a)–(8c), respectively.
B¯ 2 ¯ 12 D
Lˆ C¯ 2
1
¯ 21 . D
ˆ σ¯ F + GLH <1
(9)
where σ¯ (·) denotes the maximum singular value of a matrix, and F , G and H can be defined from the problem data and computed quantities as 1
1
1
1
γ − 2 X 2 B¯ 1 F , 1 1 γ − 2 C¯ 1 X − 2 γ −1 D11 1 1 ¯ 21 . H , C¯ 2 X − 2 γ − 2 D ¯ −2 X 2 AX
G,
1
X 2 B¯ 2
1 γ − 2 D¯ 12
(10)
For notational convenience, it is denoted that F ∈ RmF ×nF and Lˆ ∈ RmL ×nL . This paper considers the most typical case in which G is a full rank ‘‘thin’’ matrix with more rows than columns. Similarly, H is assumed to be a full rank ‘‘fat’’ matrix. The rest of this subsection reduces (9) into a simplified matrix approximation constraint to be defined in (13). First denote the SVD of G as
NG
SG 0
VG T = UG SG VG T
(11)
where UG , VG and NG are the standard orthonormal basis matrices for the column space, row space and left null space of G respectively. SG is a positive diagonal matrix of the singular values of G. Similarly, the SVD of H is
T
Bcl . Dcl
+
with blkdiag X − 2 , γ − 2 I . Then (8b) becomes
G = UG
B¯ 1 D11
Substitute this expression into (8b) and multiply both sides of (8b)
A¯ C¯ 1
right-hand-side in (8b) becomes
Rank minimization and order minimization are not the same because the inequality above can be strict. However, for any controller Kˆ there exists a controller matrix Lˆ whose dimension is the same as the full matrix L but the rank is order(Kˆ ) + min{nu , ny }.
1071
H = UH
SH 0
VH
NH
T
= UH SH VH T .
(12)
Then (Sou & Rantzer, 2010c) states that if (9) is feasible, which is the case in this paper, then the following two positive definite matrices can be defined: ∆G ,
InF − F T NG NG T F
−1
and ∆H ,
1072
K.C. Sou, A. Rantzer / Automatica 48 (2012) 1069–1076
ImF − FNH NH T F T
−1
. In addition, Sou and Rantzer (2010c) states
that (9) is equivalent to the following inequality
σ¯ Fˇ + Lˇ < 1
(13)
where Fˇ , UG T ∆H UG
− 21
UG T ∆H FVH VH T ∆G VH
12
Lˇ , (P1 )−1 Lˆ (P2 )−1 P1 , VG (SG )−1 UG T ∆H UG
P2 , VH T ∆G VH
− 12
− 21
(SH )−1 UH T
(14)
with data matrices computed in (11) and (12).
5.1. Characterizing the optimal X In Section 4 it is shown that the minimum rank of Lˆ in (8a) and (8b) is sve(Fˇ ) (cf. (16) and the discussion which follows). It turns out, from Sou and Rantzer (2010c) Theorem 1, that sve(Fˇ ) = sve(F ), where F is defined in (10) as part of the data of the problem in (15). Since F is a function of X , it is then desirable to optimize (i.e. minimize) sve(F ) with respect to X subject to (8c). This will be the main issue for the rest of this subsection. Using Sylvester’s law of inertia (e.g. Horn and Johnson (1990, p. 223)), sve the number of nonnegative eigen (F ) can be seen as 1
1
1
blkdiag X 2 , γ 2 I
1
(F T F − I ) blkdiag X 2 , γ 2 I . Alternatively, for any given X , sve(F ) is the minimum objective values of
value of the following problem. 4.2. SVD solution for generalized matrix approximation
sve (F ) = min {rank(Y )}
with
Y∈ Y
Using the equivalence between (9), (13) and (8b), as well as the fact that the matrices P1 and P2 in (14) are invertible, the proposed problem in (8a) and (8b) can be shown to be equivalent to the following two problems (i.e. having one-to-one correspondent optimal solutions with the same rank). The generalized matrix approximation problem
minimize rank Lˆ Lˆ ˆ subject to σ¯ F + GLH <1
(15)
and the classical matrix approximation problem
X
×
0
(16)
subject to σ¯ Fˇ + Lˇ < 1 where Fˇ is defined in (14). The problem in (16) can be solved using SVD via the theorem by Eckart–Young–Schmidt–Mirsky (e.g. Horn and Johnson (1990, pp. 427–428)). Define sve(Fˇ ), the singular value excess of Fˇ , as the number of singular values of Fˇ which are greater than or equal to one. Then the minimum rank of (15) (as well as (8a) and (8b)) is sve(Fˇ ). In addition, a minimum rank controller
sve(Fˇ )
matrix Lˆ can be constructed as Lˆ = −P1 σi (Fˇ )ui vi T P2 , i =1 where ui and vi are, respectively, the left and right singular vectors associated with the ordered singular values σi (Fˇ ), and P1 and P2 are defined in (14). Also, since the optimal solution to (16) is not unique, it is possible to optimize over the set of all minimum rank ˆ for instance, to minimize the controller feedthrough term Dˆ k . L, Note that even if the rank/dimension assumption on G and H is not satisfied, the results in this section still apply with a slight ˜ G˜ T where modification. For example, suppose by SVD G = GV
˜ ∈ Rm×q is full rank with m > q, and VG˜ ∈ RmX ×q is full rank G with mX > q and VG˜ T VG˜ = I. Then a variant of (15) with G replaced
˜ should be solved. Let L˜ ⋆ be an optimal solution to the variant, by G then Lˆ ⋆ = VG˜ L˜ ⋆ is optimal to (15). 5. Optimal KYP matrix X This section returns to the issue of finding an X satisfying
(8c), so that the rank minimization problem in (8a) and (8b) can be set up. (8c) can be solved as a LMI feasibility problem (e.g. using SeDuMi/Yalmip (Sturm, 1999; Löfberg, 2004) or LMILAB). Alternatively, (8c) can be written as a Riccati inequality and two LMI’s (Zhou et al., 1995). However, it turns out that the ‘‘optimal’’ (to be explained shortly) X can be characterized exactly, provided that an associated non-convex optimization problem is solved.
B¯ 1 D11
0 ¯ 1 A I C¯ 1
γ
T
B¯ 1 X − 0 D11
0 . γI
(17)
Note that in (17) if X , in addition to Y , is treated as a decision variable, then the minimization problem of sve(F ) with respect to X can be obtained. X ,Y
Lˇ
A¯ ¯C1
Y ≽ 0 s.t. Y ≻
minimize
minimize rank Lˇ
Y,
subject to
rank(Y ) X satisfies (8c) Y ∈ Y in (17), parameterized by X .
(18)
In a common heuristics to solve (18), rank(Y ) in the objective is replaced with trace(Y ) (e.g. Fazel, Hindi, and Boyd (2004)), resulting in a tractable semidefinite program (Boyd, El Ghaoui, Feron, & Balakrishnan, 1994). This is however found to be insufficient, as demonstrated in Section 6. Therefore, further improvement over the heuristic candidate X is desirable. 5.2. Updating X with a given initial guess It is possible to update X from an initial guess (e.g. the one provided by the heuristics in Fazel et al. (2004)) in such a way that the next iterate also satisfies (8c), and additionally the corresponding sve(F ) does not increase (and possibly decreases). The detail is as follows. Suppose X0 satisfies (8c). Then define R0 ∈ R(n+nk +nw )×(n+nk +nw ) such that A¯ R0 , ¯ C1
B¯ 1
D11
T
X0 0
0 ¯ 1 A I C¯ 1
γ
B¯ 1 X − 0 0 D11
0 . γI
(19)
Let of R0 in (19) be R0 = the eigenvalue decomposition T T i∈I+ λi vi vi + j∈I− λj vj vj , where λi ≥ 0 for i ∈ I+ are the nonnegative eigenvalues and λj < 0 for j ∈ I− are the negative eigenvalues. By (17), the number of elements in the index set I+ is sve(F ) for the F (X0 ). The basic idea of the local update is to ‘‘let the nonnegative eigenvalues become negative, at the expense of possibly increasing the values of the positive eigenvalues’’. For a ̸=⋆ concrete example, let i⋆ = argmini ∈I+ {λi } and I+ be such that ̸=⋆
{i⋆ } and I+ form a partition of I+ . That is, R0 = λi⋆ vi⋆ vi⋆ T + λi vi vi T + λj vj vj T . ̸=⋆ i ∈I+
j ∈I−
(20)
K.C. Sou, A. Rantzer / Automatica 48 (2012) 1069–1076
Then a particular local update problem can be formulated as the following semidefinite program with decision variables being X and c ∈ Rn+nk +nw . min ci⋆ X ,c
A¯ s.t. ¯ C1
T
B¯ 1
D11
X
≼ ci⋆ vi⋆ vi⋆ + T
0 ¯ 1 A I C¯ 1
0
γ
ci vi vi T +
̸ ⋆ = i ∈I+
X B¯ 1 − 0 D11
0 γI
cj vj vj T
j ∈I−
cj < 0 ∀j ∈ I−
X 0
0 γI
≻
Acl Ccl
Bcl Dcl
T
X 0
0 1 Acl Ccl I
γ
Bcl . Dcl
(21)
The optimal solution to (21) satisfies the following properties, whose proof is in the Appendix. Proposition 3. Let (21) be defined with X0 (cf. (19) and (20)). If X˜ is an optimal solution to (21), then it satisfies (8c) (as one of the constraints in (21)). In addition, it holds that sve(F (X˜ )) ≤ sve(F (X0 )). Therefore, successively applying the local update such as (21) results in a sequence of admissible X with non-increasing sve(F (X )). 6. Numerical assessment of the proposed reduction method In this section detailed numerical studies are reported. This is to complement the partial theoretical results concerning the restrictions introduced in Sections 3.2 and 3.3. Section 6.1 provides an example demonstrating that the proposed method has the potential to be as good as weighted balanced truncation in the applications where the later method traditionally dominates. Section 6.2 compares the proposed method against common controller reduction methods in some well-known benchmarks. This study demonstrates the advantages and drawbacks of the proposed method, and provides some guidelines of when the proposed method might work well. Section 6.3 discusses the importance of choosing the appropriate ‘‘KYP matrix’’ X . The corresponding numerical study demonstrates the importance of the heuristic update procedure described in Section 5.2. 6.1. Reduction of a Youla parameterized controller In this example the feedback setup in Fig. 1 is considered. A stable third order open loop plant P is taken from Garpinger (2009, p. 21 and 77). The plant P has 3 inputs and 3 outputs (and the controller is SISO). In Garpinger (2009), P is converted into discrete-time and then a stable controller is designed using the Youla optimized method in Boyd, Barratt, and Norman (1990), Wernrud (2008). The controller aims to, among satisfying other constraints, minimize the (time-domain) integral of the absolute value of the first closed loop output (performance measurement) when the first closed loop input (disturbance load) is a unit step function (see Garpinger (2009, p. 18) for more detail). In this example, the objective of the reduced controller is to mimic the closed loop input/output relationship of the full controller as much as possible. Therefore, the proposed controller reduction method should be applied to the feedback setup of a controller K with an ¯ which is defined such that for any reduced augmented plant P, ˆ controller K , it holds that P¯ ⋆ Kˆ , P ⋆ K − P ⋆ Kˆ .
1073
The controller has 152 states and it includes an integrator, which is customary to keep in the reduced controller. However, since the closed loop low frequency error would still be very small (cf. Fig. 2) even if the controller integrator is replaced with some poles very closed to the unit circle, the proposed method is not guaranteed to keep the controller integrator. To ensure that the reduced controller also has an integrator, the controller is (additively) divided into K = K˜ + Kint where Kint contains the integrator pole. Then Kint is added to P¯ as a feedback and only the K˜ part of the controller is reduced. After that, Kint is added back to the intermediately reduced controller to form the final one with an integrator. To apply the proposed controller reduction scheme in (8a)–(8c). The bounded-real matrix X is computed by solving a Riccati equation as described in the beginning of Section 5. This is because it is the only affordable approach to compute X considering its dimension. Applying the proposed reduction method results in a reduced controller with 18 states. As a comparison, two other reduced controllers of order 18 are obtained. One is by directly applying balanced truncation (BT) (see, e.g. Obinata and Anderson (2001, pp. 7–8)) to the full controller K , and the other is by frequency-weighted balanced truncation (FWBT) (Obinata & Anderson, 2001, pp. 106–107) with the input and output weights being (I − P0 (2, 2)K )−1 P0 (2, 1) and (I − P0 (2, 2)K )−1 P0 (1, 2) respectively (to minimize first order closed loop difference). Fig. 2 shows the magnitude Bode diagram of the error of closed loop systems due to different reduced controllers, as well as the magnitude Bode diagram of the full closed loop system for reference. The left and right figures correspond to the load sensitivity function (from input 1 to output 1) and sensitivity function (from input 1 to output 2) respectively (see Garpinger (2009, pp. 21–22) and Åström and Murray (2009, pp. 316–317) for detail). It is observed that the approximation qualities of the three controller reduction methods are similar in this example. Finally, the integrals of the absolute value of output 1 (i.e. the objective the Youla controller seeks to minimize) due to the full controller, proposed reduced controller and the balanced truncated controllers are 0.9687, 0.9686, 0.9753 and 0.9702 respectively. The lower objective value of the proposed controller is probably due to violation of some design constraint or the inexact (FIR based (Wernrud, 2008)) optimization of the Youla controller. 6.2. Reduction instances from public sources In this subsection some previously studied instances of the controller reduction problem in the setup of Fig. 1 are considered. The data of the plants P are available, in continuous-time state space matrices, from known sources. There are four instances considered, including HIMAT (an aircraft model, data available from the MathWorks website), ROC1, ROC2 and ROC3 from Complib (Leibfritz, 2004). In the experiment, three controller reduction schemes are evaluated. They include the proposed SVD based method, HIFOO (Gumussoy, Henrion, Millstone, & Overton, 2009) and the MATLAB’s balanced truncation methods (taking the best result of ncfmr and balred). ncfmr is an implementation of normalized coprime factors balanced truncation (NCFBT). On the other hand, balred is the basic balanced truncation (BT). Other (weighted) balanced truncation schemes have been reported (Zhou et al., 1995), but only ncfmr and balred are readily available in MATLAB. For the HIFOO test, HIFOO directly designs reduced order controllers from the plants P, without using the information of any full controller. HIFOO is applied five times (with the previous output being the next initial guess) for the HIMAT, ROC1 and ROC2
1074
K.C. Sou, A. Rantzer / Automatica 48 (2012) 1069–1076
(a) Load sensitivity function error magnitude.
(b) Sensitivity function error magnitude.
Fig. 2. Errors of the sensitivity functions in the example in Section 6.1. The sensitivity functions due to the full controller are also shown for reference. Table 1 Closed loop H∞ norm for different controller reduction schemes in the ROC2 example. The full order controller has order 9, with a closed loop H∞ norm of 0.04104. Method
SVD
HIFOO
NCFBT/BT
Order = 4 Order = 6 Order = 7
0.04209 0.04121 0.04104
0.04574 0.04386 0.04240
0.06141 0.04396 0.04105
Table 2 Computational time for different controller reduction schemes in the ROC2 example.
Table 4 Computational time for different controller reduction schemes in the ROC3 example. Method
SVD
HIFOO (s)
NCFBT/BT (s)
Order = 1 Order = 6 Order = 7
N/A 166 s 110 s
702 2921 3383
≈0 ≈0 ≈0
Table 5 Closed loop H∞ norm for different controller reduction schemes in the HIMAT example. The full order controller has order 10, with a closed loop H∞ norm of 0.7885.
Method
SVD (s)
HIFOO (s)
NCFBT/BT (s)
Method
SVD
HIFOO
NCFBT
Order = 4 Order = 6 Order = 7
112 151 144
2043 2964 2800
≈0 ≈0 ≈0
Order = 3 Order = 8
N/A 0.8242
1.045 0.8487
Unstable Unstable
Table 3 Closed loop H∞ norm for different controller reduction schemes in the ROC3 example. The full order controller has order 8, with a closed loop H∞ norm of 45.70. Method
SVD
HIFOO
NCFBT/BT
Order = 1 Order = 6 Order = 7
N/A 52.88 45.71
Unstable 8744 1920
Unstable Unstable 46.78
instances, and ten times for the ROC3 instance. Iterative application of HIFOO is suggested (e.g. Henrion (2006)) because HIFOO employs randomized local searches. For the NCFBT/BT test, full controllers are designed as optimal H∞ norm controllers. Then NCFBT/BT is used to construct reduced controllers. To test the proposed scheme, the plants P are first discretized with a sampling time of 0.01 s (using c2d with the tustin option in MATLAB), and then full order discrete-time optimal controllers are designed. Then the proposed scheme generates discrete-time reduced controllers, which are converted back to continuous-time to evaluate the closed loop H∞ norm. The numerical experiments exhibit advantages, as well as limitations, of the proposed controller method when compared with HIFOO and NCFBT/BT. The ROC2 and ROC3 test cases show the potential of the proposed method. Tables 1 and 2 summarize the closed loop H∞ norm performance and CPU time for the ROC2 case. Tables 3 and 4 analogously summarize the ROC3 case. When compared with HIFOO in these examples, the proposed scheme produces more accurate reduced controllers in a much shorter time (the contrast in the ROC3 case is obvious). When compared with NCFBT/BT, the proposed method generates more
Table 6 Closed loop H∞ norm for different controller reduction schemes in the ROC1 example. The full order controller has order 8, with a closed loop H∞ norm of 1.127. Method Order Order Order Order
=1 =2 =6 =7
SVD
HIFOO
NCFBT
N/A N/A 1.281 1.129
1.426 1.244 1.113 1.196
4.83 35.91 Unstable Unstable
accurate reduced controllers, but balanced truncation requires less time to compute. However, the proposed method fails to find any first order reduced controller for the ROC3 case. This is not surprising though, because the proposed method is a local search in nature. The magnitude Bode diagrams of the closed loop system in the ROC2 example (with reduced controllers of order 4) are shown in Fig. 3. These Bode diagrams depict some sensitivity functions of an aircraft model described in Gangsaas, Bruce, Blight, and Ly (1986). It can be seen that the sensitivity functions due to the proposed reduced controller are the best low order substitutes for the full order ones in this example. On the other hand, the HIMAT and ROC1 cases (in Tables 5 and 6) display the limitation of the proposed method. For example, in the HIMAT case with the current implementation, the proposed method cannot generate reduced controllers of order lower than 8 (the full order being 10). On the contrary, HIFOO manages to find a third order reduced controller with closed loop H∞ norm being 1.0451. For certain applications this may be satisfactory. In addition, in the ROC1 case HIFOO finds acceptable reduced controllers of order one and two, while the proposed
1075
Magnitude
Magnitude
K.C. Sou, A. Rantzer / Automatica 48 (2012) 1069–1076
Frequency (rad/s)
Frequency (rad/s)
(b) Sensitivity function magnitude from vertical gust noise input wn to normal vertical acceleration nz .
Magnitude
Magnitude
(a) Sensitivity function magnitude from horizontal gust input un to normal vertical acceleration nz .
Frequency (rad/s)
Frequency (rad/s)
(c) Sensitivity function magnitude from servo noise input δn to normal vertical acceleration nz .
(d) Sensitivity function magnitude from pitch rate sensor noise qn to normal vertical acceleration nz .
Fig. 3. Sensitivity functions due to the full controller and reduced controllers of order four in the ROC2 example.
method cannot find any. On the other hand, balanced truncation deliver quite lackluster performance in these test cases. In summary, the examples show that the proposed controller reduction scheme manages to find the middle ground between a global optimization based approach (i.e. HIFOO) and local approximation approach (i.e. NCFBT/BT). Hence, it can be useful in complementing the existing schemes. 6.3. Further discussions on algorithmic details Finally, note that the proposed method is applied in Section 6.2 iteratively. First, with the given full controller matrix L, the ‘‘trace heuristic’’ version of problem (18) in Section 5.1 is solved for an initial guess of X . Then, the local update problem in (21) in Section 5.2 is solved several times (e.g. three in this example) for further reduction of the sve(F (X )) (pertaining to the order of the reduced controller). After X is determined, the SVD based ˆ The above process, computation in Section 4 is applied to find L. ˆ from L to L, can be repeated until no more significant progress can be made. The computation time for the SVD method in Tables 2 and 4 include all of the operations above, for all the iterations. It is instructive to see the effect of the local update (21) in Section 5.2, by comparing the achievable reduced controller order with and without the update procedure. Without using the update procedure, for the ROC2 example the reduced controllers generated by the proposed method would have orders of 9, 9, and 8 (as opposed to 4, 6 and 7 in Table 1). For the ROC3 example, the corresponding reduced orders would be 7 and 8 (as opposed to 6 and 7 in Table 3). This demonstrates that computation effort should be spent on solving (18) to as close to optimality as possible.
7. Conclusion This paper demonstrates that it is possible to obtain a closed loop stable controller reduction scheme which linearizes the bounded-real lemma using full controller information—in the spirit of balanced truncation. The proposed controller reduction problem is a minimum rank generalized matrix approximation problem, which can be solved efficiently using SVD. In this paper, numerical examples have demonstrated that proposed controller reduction scheme is on par with the more well-known alternatives such as balanced truncation and HIFOO. Finally, a similar controller reduction setup in H2 norm is also possible. See Sou and Rantzer (2010b) for detail. Acknowledgments Support from the Swedish Research Council through the Linnaeus Centers ACCESS and LCCC are gratefully acknowledged. The authors would also like to thank the anonymous reviewers for helpful comments. Appendix A.1. Proof of Proposition 2 First prove the inequality order(Kˆ )
≤ rank Aˆ k Bˆ k . Bˆ k . The inequality
ˆk , rank A ˆ order(K ) ≤ r holds trivially when r ≥ nˆ k . Therefore, for the rest For convenience, denote r
1076
K.C. Sou, A. Rantzer / Automatica 48 (2012) 1069–1076
of the proof it is assumed that r < nˆ k . Because r < nˆ k , there exists (e.g. by SVD) a unitary matrix S ∈ Rnˆ k ׈nk such that S T S = Inˆ k and S T Aˆ k
Bˆ k = S T Aˆ k
S T Bˆ k =
A˜ k1 0A
B˜ k1 0B
where the zero matrices 0A and 0B have at least nˆ k − r rows. Thus, Kˆ defined by Lˆ is a similarity transform, with transformation matrix S, of a LTI system K˜ , whose system matrices are
A˜ k1 S 0A
,
B˜ k1 0B
, Cˆ k S
ˆ k respectively. Since K˜ has at least nˆ k − r uncontrollable and D states, order(K˜ ) ≤ r. Finally, since Kˆ and K˜ have the same order as they are similarity transforms, it is concluded that order(Kˆ ) ≤ r = rank Aˆ k Bˆ k . The proof of the inequality order(Kˆ ) ≤ rank
ˆ Ak Cˆ k
is omitted because it is similar to the
previous case. Finally, combining the two inequalities leads to the desired result. A.2. Proof of Proposition 3 To show that sve(F (X˜ )) ≤ sve(F (X0 )), denote, for any symmetric matrix M , nne(M ) as the number of nonnegative eigenvalues of M. Then (17) and the discussion immediately before it implies that sve(F (X˜ )) is equal to
A¯ nne ¯ C1
B¯ 1 D11
T
X˜
0
0 ¯ A 1 I C¯ 1
γ
˜ B¯ 1 − X 0 D11
0 γI
(22)
(22) is less than or equal to
nne ci⋆ vi⋆ vi⋆ T +
ci vi vi T +
̸ ⋆ = i ∈I+
j ∈I−
cj vj vj T .
(23)
This relation can be inferred by applying the first constraint in (21) to (17). (23) is less than or equal to
nne λi⋆ vi⋆ vi⋆ T +
̸ ⋆ = i ∈I+
λ i vi vi T +
j ∈I−
λj vj vj T
(24)
and the inequality is true because of the following. X0 , together with the λi ’s (i.e. the eigenvalues of R0 in (20)), is feasible for (21). ̸=⋆ Therefore, for all j ∈ I− , both cj and λj are negative. For all i ∈ I+ , both ci and λi are nonnegative. However, it might be possible that ci⋆ < 0 ≤ λi⋆ , and hence this inequality can be strict. Finally, (24) is equal to sve(F (X0 )) because of (19), (20) and (22). This concludes that sve(F (X˜ )) ≤ sve(F (X0 )). References Ajala, O., Schuler, S., & Allgower, F. (2010). Linf-Gain controller order reduction for discrete-time systems. In American control conference. Åström, K., & Murray, R. (2009). Feedback systems: an introduction for scientists and engineers. Princeton University Press. Anderson, B. D. O., & Liu, Y. (1989). Controller reduction: concepts and approaches. IEEE Transactions on Automatic Control, 34(8). Ankelhed, D., Helmersson, A., & Hansson, A. (2006). Suboptimal model reduction using LMIs with convex constraints, Technical report. Department of Electrical Engineering. Linköping University. Boyd, S., Barratt, C., & Norman, S. (1990). Linear controller design: limits of performance via convex optimization. Proceedings of the IEEE, 78(3), 529–574. Boyd, S., El Ghaoui, L., Feron, E., & Balakrishnan, V. (1994). Linear matrix inequalities in system and control theory. SIAM. Doyle, J., Glover, K., Khargonekar, P., & Francis, B. (1989). State-space solutions to standard H2 and H∞ control problems. IEEE Transactions on Automatic control, 34(8), 831–847.
Enns, D. (1984). Model reduction with balanced realizations: an error bound and a frequency weighted generalization. Conference on Decision and Control. Fazel, M., Hindi, H., & Boyd, S. (2004). Rank minimization and applications in system theory. In American control conference. Boston. Gahinet, P., & Apkarian, P. (1994). A linear matrix inequality approach to H∞ control. International Journal of Robust and Nonlinear Control, 4, 421–448. Gangsaas, D., Bruce, K., Blight, J., & Ly, Uy-Loi (1986). Application of modem synthesis to aircraft control: three case studies. Automatic Control, IEEE Transactions on, 31(11), 995–1014. Garpinger, O. (2009). Design of robust PID controllers with constrained control signal activity. Licentiate Thesis. Reglerteknik. Lund University. Sweden. Goddard, P.J., & Glover, K. (1993). Controller reduction: weights for stability and performance preservation. In Proccedings of 32nd IEEE conference on decision and control (pp. 2903–2908). Gumussoy, S., Henrion, D., Millstone, M., & Overton, M.L. (2009). Multiobjective robust control with HIFOO 2.0. In Proceedings of the IFAC symposium on robust control design. Haifa, Israel. Henrion, D. (2006). Some control design experiments with HIFOO. Technical report 06570. LASS-CNRS. Horn, R., & Johnson, C. (1990). Matrix analysis. Cambridge University Press. Iwasaki, T., & Skelton, R. E. (1994). All controllers for the general H∞ control problem: LMI existence conditions and state space formulas. Automatica, 30(8), 1307–1317. Iwasaki, T., & Skelton, R. E. (1995a). The XY -centring algorithm for the dual LMI problem: a new approach to fixed-order control design. International Journal of Control, 62(6), 1257–1272. Iwasaki, T., & Skelton, R. E. (1995b). A unified approach to fixed-order controller design via linear matrix inequalities. Mathematical Problems in Engineering, 1(1), 59–75. Leibfritz, F. (2004). COMPleib: COnstraint matrix-optimization problem library—a collection of test examples for nonlinear semidefinite programs. Control system design and related problems. Technical report. University of Trier. Department of Mathematics. data available from http://www.complib.de/. Obinata, G., & Anderson, B. D. O. (2001). Model reduction for control system design. Springer. Sandberg, H. (2010). An extension to balanced truncation with application to structured model reduction. IEEE Transactions on Automatic Control, 55(4), 1038–1043. Sou, K.C., & Rantzer, A. (2010a). A singular value decomposition based closed loop stability preserving controller reduction method. American Control Conference, Baltimore, 2010, pp. 1079–1084. Sou, K.C., & Rantzer, A. (2010b). Controller reduction with closed loop error guarantee. 49th Conference on Decision and Control, Atlanta, 2010. Sou, K.C., & Rantzer, A. (2010c). On the minimum rank of a generalized matrix approximation problem in the maximum singular value norm. In 19th international symposium on mathematical theory of networks and systems. Sturm, J. F. (1999). Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones. Optimization Methods and Software, 11–12, 625–653. Wernrud, A. (2008). QTool 0.1 — Reference manual. Department of automatic control. Lund University, Sweden. Löfberg, J. (2004). Yalmip: a toolbox for modeling and optimization in MATLAB. In Proceedings of the CACSD Conference. Taipei, Taiwan. Zhou, K. (1993). Frequency weighted model reduction with L∞ error bounds. Systems and Control Letters, 21, 115–125. Zhou, K., Doyle, J., & Glover, K. (1995). Robust and optimal control. Prentice Hall.
Kin Cheong Sou received a Ph.D. degree in Electrical Engineering and Computer Science at Massachusetts Institute of Technology in 2008. From 2008 to 2010 he was a postdoctoral researcher at Lund University, Lund, Sweden, and since September 2010 he has been a postdoctoral researcher at KTH Royal Institute of Technology, Stockholm, Sweden. His research interests include power system cybersecurity analysis, environment aware building and community, convex/non-convex optimization and model reduction for dynamical systems. Anders Rantzer received a Ph.D. degree in Optimization and Systems Theory in 1991 from the Royal Institute of Technology (KTH), Stockholm, Sweden. After postdoctoral positions at KTH and at IMA, University of Minnesota, he joined the Department of Automatic Control at Lund University in 1993. He was appointed professor of Automatic Control in Lund 1999. The academic year of 2004/05 he was a visiting associate faculty member at Caltech. Since 2008 he has coordinated the Linnaeus center LCCC at Lund University. Rantzer has been serving as an associate editor of IEEE Transactions on Automatic Control and several other journals. He is a winner of the SIAM Student Paper Competition, the IFAC Congress Young Author Price and the IET Premium Award for the best article in IEE Proceedings — Control Theory & Applications during 2006. He is a Fellow of IEEE and a member of the Royal Swedish Academy of Engineering Sciences. His research interests are in modeling, analysis and synthesis of control systems, with particular attention to uncertainty, optimization and distributed control.
Contents lists available at SciVerse ScienceDirect
Automatica journal homepage: www.elsevier.com/locate/automatica
Controller reduction via minimum rank matrix approximation✩ Kin Cheong Sou a,1 , Anders Rantzer b a
Automatic Control Lab and ACCESS Linnaeus Center, School of Electrical Engineering, Royal Institute of Technology, Osquldas väg 10, floor 6, SE-100 44 Stockholm, Sweden
b
Department of Automatic Control, Lund University, Box 118, SE-221 00 Lund, Sweden
article
abstract
info
Article history: Received 7 April 2011 Received in revised form 17 August 2011 Accepted 27 October 2011 Available online 30 March 2012
In this paper a controller reduction method for discrete-time linear time-invariant systems is described. Using the bounded-real lemma, the proposed method generates reduced controllers with closed loop stability and H∞ norm performance guarantee. Information of the full controller is used as a basis for reduction using singular value decomposition. This is different from traditional model reduction schemes such as weighted balanced truncation. Numerical assessment of the proposed method is given in the end of the paper. © 2012 Elsevier Ltd. All rights reserved.
Keywords: Model reduction Singular value decomposition
1. Introduction In this paper a controller reduction problem for discrete-time linear time-invariant systems depicted in Fig. 1 is considered. In this reduction problem, the feedback structure should remain and only the controller can be simplified. This can be viewed as a special case of the more general problem of structured model reduction (e.g. Sandberg (2010)). However, the controller reduction problem, by itself, has also been extensively studied because of its potential values. The main idea of the proposed controller reduction method can be explained as follows. Standard H∞ and H2 synthesis procedures involve solution of the following linear matrix inequality (LMI):
X 0
0 γI
A ≻ cl Ccl
Bcl Dcl
T
X
0
0 1 Acl Ccl I
γ
Bcl Dcl
(1)
where Acl , Bcl , Ccl and Dcl are affine functions of an unknown matrix L describing the controller, and X is another unknown as the ‘‘KYP certificate’’. The optimal L and X found by solving (1) often lead to high order controllers. We are interested in finding good approximations of low order. The approach is to use the same X as
✩ The material in this paper was partially presented at the 2010 American Control Conference, June 30–July 2, 2010, Baltimore, Maryland, USA and the 49th IEEE Conference on Decision and Control (CDC), December 15–17, 2010, Atlanta, Georgia, USA. This paper was recommended for publication in revised form by Associate Editor Richard D. Braatz under the direction of Editor Frank Allgöwer. E-mail addresses: [email protected] (K.C. Sou), [email protected] (A. Rantzer). 1 Tel.: +46 8 7907427; fax: +46 8 7907329.
0005-1098/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.automatica.2012.03.011
for the optimal high order controller, then replace L by a matrix Lˆ of minimal rank, subject to the inequality (1). Rank minimization is not equivalent to controller order reduction, but low rank leads to low order controllers and works well in examples. The above idea can be viewed as the middle ground between two existing groups of controller reduction methods, balancing computation effort and approximation quality. One group is balanced truncation/Hankel norm reduction and their weighted variants (e.g. Enns (1984), Anderson and Liu (1989), Goddard and Glover (1993), Zhou (1993), and Obinata and Anderson (2001)). These methods utilize full controller information. They are efficient and produce good local approximation results. The other group of known reduction methods is optimization based, typically formulated using the KYP lemma (e.g. Gahinet and Apkarian (1994), Iwasaki and Skelton (1994), Ankelhed, Helmersson, and Hansson (2006), Ajala, Schuler, and Allgower (2010), and Iwasaki and Skelton (1995a,b)). These methods involve more expensive computations, but a globally optimal reduced controller can be found if the corresponding optimization problems are solved to optimality. Full controller information is not always necessary in these optimization based methods. The purpose of this paper is not to present yet another attempt to solve the KYP problem to optimality. Rather, it describes a perspective on how to utilize the full controller information to obtain a reduced controller. This perspective is different from weighted balanced truncation, and it attempts to combine the benefits of weighted balanced truncation and KYP optimization methods. In addition, via numerical examples this paper provides some assessment of how well the investigated perspective can work in practice. Nevertheless, the full theoretical characterization of the limit of performance remains unknown in this paper.
1070
K.C. Sou, A. Rantzer / Automatica 48 (2012) 1069–1076 w
w
z P
u
y
z P
u
K
With the matrices defined in (4), the state space matrices of the reduced closed loop system P ⋆ Kˆ can be described as affine functions of the reduced controller Lˆ as
y
ˆ cl , A¯ + B¯ 2 Lˆ C¯ 2 A
K
Fig. 1. The block diagram on the left is the feedback interconnection of a plant P and a controller K considered in this paper. The problem of controller reduction is to find a lower order controller Kˆ to replace K , so that the performances of the reduced closed loop system (on the right) do not degrade too much.
The rest of the paper is organized as follows. In Section 2 the background material and notations are introduced. In Section 3 the proposed controller reduction problem is motivated in detail and then formulated. Then Section 4 describes the SVD based solution to the proposed problem. The order of the reduced controller turns out to be dependent on a ‘‘KYP matrix’’, the search of which is to be discussed in Section 5. Numerical assessments of the proposed method are presented in Section 6. 2. Background and notation 2.1. Controller reduction problem In this paper, all the systems involved (except some in Section 6) are discrete-time linear time-invariant systems. The main object of interest is the feedback setup of a plant P and a controller K in Fig. 1. The feedback closed loop system is denoted as P ⋆ K , and it is assumed to be stable. In the controller reduction problem in this paper, P and K are given data. A reduced controller Kˆ is sought as the minimum order controller whose corresponding closed loop system, denoted as P ⋆ Kˆ , has H∞ norm less than a userspecified upper bound. The controller reduction problem as stated is computationally difficult to solve, and therefore an approximate solution will be found by solving a different but related problem proposed in Section 3. 2.2. State space characterization of closed loop system In Fig. 1 the dimensions of the signals w, u, z , y are nw , nu , nz , ny respectively. The plant P is described in state space matrices as follows.
C1 P = (zI − A)−1 B1 C2
B2
D + 11 D21
D12 D22
(2)
where A ∈ Rn×n , D12 ∈ Rnz ×nu , D21 ∈ Rny ×nw and z ∈ C. In this paper, a technical assumption is made that D22 = 0 (see Doyle, Glover, Khargonekar, and Francis (1989) for detail). The given full controller K is described by its state space matrices as K = Ck (zI − Ak )−1 Bk + Dk . Similarly, the reduced ˆ k . In above, Ak ∈ Rnk ×nk , Dk ∈ controller Kˆ is Cˆ k (zI − Aˆ k )−1 Bˆ k + D ˆ nu ×ny n ׈ n n × ˆ k ∈ R u ny . The subsequent exposition R and Aˆ k ∈ R k k , D requires the definitions of the following two matrices.
Ak L, Ck
Bk Dk
Aˆ k ˆ and L , Cˆ k
Bˆ k ˆk . D
(3)
Subsequently, L and Lˆ are referred to as the full and reduced controller matrices, respectively. Also, define the data matrices (A¯ , B¯ 1 , B¯ 2 , etc.) A 0 B1 0 B2 0 0 0 Inˆ k 0 B¯ 1 B¯ 2 A¯ nˆ k ¯ ¯ 12 , C1 C1 D11 D 0 D11 0 D12 . ¯ 21 ¯ 22 C¯ 2 D D 0 Inˆ 0 0 0 C2
k
0
D21
0
0
(4)
¯ 12 Lˆ C¯ 2 Cˆ cl , C¯ 1 + D
ˆ cl , B¯ 1 + B¯ 2 Lˆ D¯ 21 B
(5)
ˆ cl , D11 + D¯ 12 Lˆ D¯ 21 . D
The dimensions of the closed loop state space matrices above depend on the dimension of Lˆ (i.e. nˆ k ). Also, let Acl , Bcl , Ccl and Dcl denote the state space matrices of the full closed loop system P ⋆ K . In the case where nˆ k = nk , which is the main focus of this paper, these matrices can also be written as in (5), with L replacing ˆ Finally, ∥P ⋆ K ∥∞ and ∥P ⋆ Kˆ ∥∞ denote the H∞ norms of the L. respective closed loop systems. 3. Basic problem formulation 3.1. Closed loop performance objectives Two performance objectives should be guaranteed by the reduced controller Kˆ . First, P ⋆ Kˆ should be stable. Second, for any given γ > ∥P ⋆ K ∥∞ , it holds that ∥P ⋆ Kˆ ∥∞ < γ (the lower bound on γ ensures that K can be one of the candidate reduced controllers). These two performance objectives can be specified using the state space matrices of P ⋆ Kˆ (which, according to (5), ˆ via the bounded-real lemma as follows. are functions of L), Theorem 1 (Zhou, Doyle, & Glover, 1995). Let γ be a positive scalar. In addition, denote the closed loop system corresponding to Lˆ (cf. (3)
ˆ cl and (5)) as P ⋆ Kˆ = Cˆ cl zI − A
−1
ˆ cl + Dˆ cl . The following two B
statements are equivalent.
• P ⋆ Kˆ is stable and ∥P ⋆ Kˆ ∥∞ < γ • there exists X ≻ 0 such that 0 ˆ ˆ cl Bˆ cl T X X 0 A 1 Acl ≻ ˆ ˆ cl 0 γI I 0 Ccl D Cˆ cl γ
(6a)
ˆ cl B ˆ cl . D
(6b)
Remark 1. The results in this paper rely on the structure in (6b). Since the continuous-time bounded-real lemma does not have a statement similar to (6b) in structure, this paper considers the discrete-time case only. 3.2. A convex restriction of controller matrix candidates It is a well-known difficulty to apply (6b) to simultaneously solve for the bound-real matrices X and reduced controller ˆ Refs. (Gahinet & Apkarian, 1994; Iwasaki & Skelton, matrices L. 1994) report a convex equivalence of (6b), but it remains nonconvex if an order constraint is imposed (i.e., nˆ k < nk ). Instead, this paper proposes to use the full controller matrix L to separate ˆ First X is found as a solution to the following the search of X and L. LMI
X 0
0 γI
A ≻ cl Ccl
Bcl Dcl
T
X
0
0 1 Acl Ccl I
γ
Bcl Dcl
(7)
where Acl , Bcl , Ccl and Dcl are defined in the discussion following (5). (7) ensures that X is admissible in the sense that at least L satisfies (6b) with X defined above. Once X is determined, the corresponding Lˆ is found by solving (6b). Imposing (7) leads to a restriction on the choice of X as opposed to letting it be an unknown in (6b). Nevertheless, often the full controller is
K.C. Sou, A. Rantzer / Automatica 48 (2012) 1069–1076
unnecessarily high order. Hence, it is sensible to expect some lower order controllers would also satisfy (6b) with X from (7) with an appropriate γ which leads to acceptable closed loop performances. In addition, imposing (7) is in line with local methods such as weighted balanced truncation — to utilize the full controller information as well as possible. In this paper, a detailed evaluation through numerical examples is conducted to assess the limit of performance imposed by (7), though a full theoretical investigation remains difficult. 3.3. Order specification via rank of the controller matrix Apart from the restrictiveness, (7) also complicates the search of controllers of reduced orders. In particular, the combination of (5) and (7) implies that nˆ k = nk (i.e. the A matrices of the full and reduced controllers have the same dimension). Hence, to achieve order reduction it would be desirable to directly ˆ However, this renders minimize order(Kˆ ) (which is a function of L). ˆ the search of L intractable. Instead, this paper proposes that rank(Lˆ ) should be minimized. This is based on the following result preliminarily from Sou and Rantzer (2010a,b) (now proved in the Appendix). Proposition 2. Let Kˆ be the reduced controller associated with a reduced controller matrix Lˆ as in (3), then
order(Kˆ ) ≤ min rank Aˆ k Aˆ k Cˆ k
≤ rank
Bˆ k ˆk D
Bˆ k
Aˆ k , rank Cˆ k
The reduced controller is found via the following steps: (1) Solve the LMI’s in (8c) to obtain X . The solutions to (8c) are typically non-unique, and there are different approaches to find the different X . The details are given in Section 5. (2) Once X is determined from (8c), solve the constrained rank minimization problem of (8a) and (8b). This problem can be solved efficiently using SVD. See Section 4. (3) Once Lˆ of the minimum rank is determined, form a minimal state space model from Lˆ as the proposed reduced order controller in the problem in Fig. 1. 4. Minimum rank controller matrix Lˆ The main issue in this section is to solve the constrained minimum rank problem (8a) and (8b), assuming that the matrix X has been determined from (8c). The main result here is that the seemingly intractable rank minimization problem in (8a) and (8b) can be shown to be equivalent to a SVD solvable classical matrix approximation problem to be defined in (16). To show the equivalence, first the matrix inequality in (8b) is simplified using Theorem 1 in Sou and Rantzer (2010c). Then, (8a) and the simplified version of (8b) becomes a classical minimum rank matrix approximation problem, solvable using SVD. 4.1. Simplification of matrix inequality The main objective of this subsection is to show that (8b) and a simplified constraint to be defined in (13) are equivalent. The development is as follows. Using on the (5), the last matrix
= rank(Lˆ ).
ˆ for the low Nevertheless, (7) might exclude the forming of such L’s order controllers. This remains an open question, but Section 6 provides an experimental indication that rank minimization need not be restrictive. 3.4. Formulation of the basic problem The data of the proposed problem are the plant matrices in (4) and L in (3) describing the full order controller. There is a user-specified scalar parameter γ such that γ > ∥P ⋆ K ∥∞ and the proposed problem would enforce that γ > ∥P ⋆ Kˆ ∥∞ , thus controlling the tradeoff between performance loss and controller order reduction. The reduced controller matrix Lˆ is the decision variable. Lˆ and L have the same dimensions (i.e. nˆ k = nk ), as discussed before. The proposed problem, in its basic form, can be formulated as minimize rank(Lˆ ) subject to
(8a)
Lˆ
X 0
0
γI
ˆ A ≻ ˆ cl Ccl
ˆ cl B ˆ cl D
T
X
0
0 ˆ 1 Acl I Cˆ cl
γ
ˆ cl B ˆ cl D
(8b)
where X ≻ 0 can be computed from the data by solving
X 0
0 γI
A ≻ cl Ccl
Bcl Dcl
T
X
0
0 1 Acl Ccl I
γ
1
(8c)
The definitions of the closed loop system matrices in (8b) and (8c) can be found in (5). Also, see Sections 3.1–3.3 for the explanations of (8a)–(8c), respectively.
B¯ 2 ¯ 12 D
Lˆ C¯ 2
1
¯ 21 . D
ˆ σ¯ F + GLH <1
(9)
where σ¯ (·) denotes the maximum singular value of a matrix, and F , G and H can be defined from the problem data and computed quantities as 1
1
1
1
γ − 2 X 2 B¯ 1 F , 1 1 γ − 2 C¯ 1 X − 2 γ −1 D11 1 1 ¯ 21 . H , C¯ 2 X − 2 γ − 2 D ¯ −2 X 2 AX
G,
1
X 2 B¯ 2
1 γ − 2 D¯ 12
(10)
For notational convenience, it is denoted that F ∈ RmF ×nF and Lˆ ∈ RmL ×nL . This paper considers the most typical case in which G is a full rank ‘‘thin’’ matrix with more rows than columns. Similarly, H is assumed to be a full rank ‘‘fat’’ matrix. The rest of this subsection reduces (9) into a simplified matrix approximation constraint to be defined in (13). First denote the SVD of G as
NG
SG 0
VG T = UG SG VG T
(11)
where UG , VG and NG are the standard orthonormal basis matrices for the column space, row space and left null space of G respectively. SG is a positive diagonal matrix of the singular values of G. Similarly, the SVD of H is
T
Bcl . Dcl
+
with blkdiag X − 2 , γ − 2 I . Then (8b) becomes
G = UG
B¯ 1 D11
Substitute this expression into (8b) and multiply both sides of (8b)
A¯ C¯ 1
right-hand-side in (8b) becomes
Rank minimization and order minimization are not the same because the inequality above can be strict. However, for any controller Kˆ there exists a controller matrix Lˆ whose dimension is the same as the full matrix L but the rank is order(Kˆ ) + min{nu , ny }.
1071
H = UH
SH 0
VH
NH
T
= UH SH VH T .
(12)
Then (Sou & Rantzer, 2010c) states that if (9) is feasible, which is the case in this paper, then the following two positive definite matrices can be defined: ∆G ,
InF − F T NG NG T F
−1
and ∆H ,
1072
K.C. Sou, A. Rantzer / Automatica 48 (2012) 1069–1076
ImF − FNH NH T F T
−1
. In addition, Sou and Rantzer (2010c) states
that (9) is equivalent to the following inequality
σ¯ Fˇ + Lˇ < 1
(13)
where Fˇ , UG T ∆H UG
− 21
UG T ∆H FVH VH T ∆G VH
12
Lˇ , (P1 )−1 Lˆ (P2 )−1 P1 , VG (SG )−1 UG T ∆H UG
P2 , VH T ∆G VH
− 12
− 21
(SH )−1 UH T
(14)
with data matrices computed in (11) and (12).
5.1. Characterizing the optimal X In Section 4 it is shown that the minimum rank of Lˆ in (8a) and (8b) is sve(Fˇ ) (cf. (16) and the discussion which follows). It turns out, from Sou and Rantzer (2010c) Theorem 1, that sve(Fˇ ) = sve(F ), where F is defined in (10) as part of the data of the problem in (15). Since F is a function of X , it is then desirable to optimize (i.e. minimize) sve(F ) with respect to X subject to (8c). This will be the main issue for the rest of this subsection. Using Sylvester’s law of inertia (e.g. Horn and Johnson (1990, p. 223)), sve the number of nonnegative eigen (F ) can be seen as 1
1
1
blkdiag X 2 , γ 2 I
1
(F T F − I ) blkdiag X 2 , γ 2 I . Alternatively, for any given X , sve(F ) is the minimum objective values of
value of the following problem. 4.2. SVD solution for generalized matrix approximation
sve (F ) = min {rank(Y )}
with
Y∈ Y
Using the equivalence between (9), (13) and (8b), as well as the fact that the matrices P1 and P2 in (14) are invertible, the proposed problem in (8a) and (8b) can be shown to be equivalent to the following two problems (i.e. having one-to-one correspondent optimal solutions with the same rank). The generalized matrix approximation problem
minimize rank Lˆ Lˆ ˆ subject to σ¯ F + GLH <1
(15)
and the classical matrix approximation problem
X
×
0
(16)
subject to σ¯ Fˇ + Lˇ < 1 where Fˇ is defined in (14). The problem in (16) can be solved using SVD via the theorem by Eckart–Young–Schmidt–Mirsky (e.g. Horn and Johnson (1990, pp. 427–428)). Define sve(Fˇ ), the singular value excess of Fˇ , as the number of singular values of Fˇ which are greater than or equal to one. Then the minimum rank of (15) (as well as (8a) and (8b)) is sve(Fˇ ). In addition, a minimum rank controller
sve(Fˇ )
matrix Lˆ can be constructed as Lˆ = −P1 σi (Fˇ )ui vi T P2 , i =1 where ui and vi are, respectively, the left and right singular vectors associated with the ordered singular values σi (Fˇ ), and P1 and P2 are defined in (14). Also, since the optimal solution to (16) is not unique, it is possible to optimize over the set of all minimum rank ˆ for instance, to minimize the controller feedthrough term Dˆ k . L, Note that even if the rank/dimension assumption on G and H is not satisfied, the results in this section still apply with a slight ˜ G˜ T where modification. For example, suppose by SVD G = GV
˜ ∈ Rm×q is full rank with m > q, and VG˜ ∈ RmX ×q is full rank G with mX > q and VG˜ T VG˜ = I. Then a variant of (15) with G replaced
˜ should be solved. Let L˜ ⋆ be an optimal solution to the variant, by G then Lˆ ⋆ = VG˜ L˜ ⋆ is optimal to (15). 5. Optimal KYP matrix X This section returns to the issue of finding an X satisfying
(8c), so that the rank minimization problem in (8a) and (8b) can be set up. (8c) can be solved as a LMI feasibility problem (e.g. using SeDuMi/Yalmip (Sturm, 1999; Löfberg, 2004) or LMILAB). Alternatively, (8c) can be written as a Riccati inequality and two LMI’s (Zhou et al., 1995). However, it turns out that the ‘‘optimal’’ (to be explained shortly) X can be characterized exactly, provided that an associated non-convex optimization problem is solved.
B¯ 1 D11
0 ¯ 1 A I C¯ 1
γ
T
B¯ 1 X − 0 D11
0 . γI
(17)
Note that in (17) if X , in addition to Y , is treated as a decision variable, then the minimization problem of sve(F ) with respect to X can be obtained. X ,Y
Lˇ
A¯ ¯C1
Y ≽ 0 s.t. Y ≻
minimize
minimize rank Lˇ
Y,
subject to
rank(Y ) X satisfies (8c) Y ∈ Y in (17), parameterized by X .
(18)
In a common heuristics to solve (18), rank(Y ) in the objective is replaced with trace(Y ) (e.g. Fazel, Hindi, and Boyd (2004)), resulting in a tractable semidefinite program (Boyd, El Ghaoui, Feron, & Balakrishnan, 1994). This is however found to be insufficient, as demonstrated in Section 6. Therefore, further improvement over the heuristic candidate X is desirable. 5.2. Updating X with a given initial guess It is possible to update X from an initial guess (e.g. the one provided by the heuristics in Fazel et al. (2004)) in such a way that the next iterate also satisfies (8c), and additionally the corresponding sve(F ) does not increase (and possibly decreases). The detail is as follows. Suppose X0 satisfies (8c). Then define R0 ∈ R(n+nk +nw )×(n+nk +nw ) such that A¯ R0 , ¯ C1
B¯ 1
D11
T
X0 0
0 ¯ 1 A I C¯ 1
γ
B¯ 1 X − 0 0 D11
0 . γI
(19)
Let of R0 in (19) be R0 = the eigenvalue decomposition T T i∈I+ λi vi vi + j∈I− λj vj vj , where λi ≥ 0 for i ∈ I+ are the nonnegative eigenvalues and λj < 0 for j ∈ I− are the negative eigenvalues. By (17), the number of elements in the index set I+ is sve(F ) for the F (X0 ). The basic idea of the local update is to ‘‘let the nonnegative eigenvalues become negative, at the expense of possibly increasing the values of the positive eigenvalues’’. For a ̸=⋆ concrete example, let i⋆ = argmini ∈I+ {λi } and I+ be such that ̸=⋆
{i⋆ } and I+ form a partition of I+ . That is, R0 = λi⋆ vi⋆ vi⋆ T + λi vi vi T + λj vj vj T . ̸=⋆ i ∈I+
j ∈I−
(20)
K.C. Sou, A. Rantzer / Automatica 48 (2012) 1069–1076
Then a particular local update problem can be formulated as the following semidefinite program with decision variables being X and c ∈ Rn+nk +nw . min ci⋆ X ,c
A¯ s.t. ¯ C1
T
B¯ 1
D11
X
≼ ci⋆ vi⋆ vi⋆ + T
0 ¯ 1 A I C¯ 1
0
γ
ci vi vi T +
̸ ⋆ = i ∈I+
X B¯ 1 − 0 D11
0 γI
cj vj vj T
j ∈I−
cj < 0 ∀j ∈ I−
X 0
0 γI
≻
Acl Ccl
Bcl Dcl
T
X 0
0 1 Acl Ccl I
γ
Bcl . Dcl
(21)
The optimal solution to (21) satisfies the following properties, whose proof is in the Appendix. Proposition 3. Let (21) be defined with X0 (cf. (19) and (20)). If X˜ is an optimal solution to (21), then it satisfies (8c) (as one of the constraints in (21)). In addition, it holds that sve(F (X˜ )) ≤ sve(F (X0 )). Therefore, successively applying the local update such as (21) results in a sequence of admissible X with non-increasing sve(F (X )). 6. Numerical assessment of the proposed reduction method In this section detailed numerical studies are reported. This is to complement the partial theoretical results concerning the restrictions introduced in Sections 3.2 and 3.3. Section 6.1 provides an example demonstrating that the proposed method has the potential to be as good as weighted balanced truncation in the applications where the later method traditionally dominates. Section 6.2 compares the proposed method against common controller reduction methods in some well-known benchmarks. This study demonstrates the advantages and drawbacks of the proposed method, and provides some guidelines of when the proposed method might work well. Section 6.3 discusses the importance of choosing the appropriate ‘‘KYP matrix’’ X . The corresponding numerical study demonstrates the importance of the heuristic update procedure described in Section 5.2. 6.1. Reduction of a Youla parameterized controller In this example the feedback setup in Fig. 1 is considered. A stable third order open loop plant P is taken from Garpinger (2009, p. 21 and 77). The plant P has 3 inputs and 3 outputs (and the controller is SISO). In Garpinger (2009), P is converted into discrete-time and then a stable controller is designed using the Youla optimized method in Boyd, Barratt, and Norman (1990), Wernrud (2008). The controller aims to, among satisfying other constraints, minimize the (time-domain) integral of the absolute value of the first closed loop output (performance measurement) when the first closed loop input (disturbance load) is a unit step function (see Garpinger (2009, p. 18) for more detail). In this example, the objective of the reduced controller is to mimic the closed loop input/output relationship of the full controller as much as possible. Therefore, the proposed controller reduction method should be applied to the feedback setup of a controller K with an ¯ which is defined such that for any reduced augmented plant P, ˆ controller K , it holds that P¯ ⋆ Kˆ , P ⋆ K − P ⋆ Kˆ .
1073
The controller has 152 states and it includes an integrator, which is customary to keep in the reduced controller. However, since the closed loop low frequency error would still be very small (cf. Fig. 2) even if the controller integrator is replaced with some poles very closed to the unit circle, the proposed method is not guaranteed to keep the controller integrator. To ensure that the reduced controller also has an integrator, the controller is (additively) divided into K = K˜ + Kint where Kint contains the integrator pole. Then Kint is added to P¯ as a feedback and only the K˜ part of the controller is reduced. After that, Kint is added back to the intermediately reduced controller to form the final one with an integrator. To apply the proposed controller reduction scheme in (8a)–(8c). The bounded-real matrix X is computed by solving a Riccati equation as described in the beginning of Section 5. This is because it is the only affordable approach to compute X considering its dimension. Applying the proposed reduction method results in a reduced controller with 18 states. As a comparison, two other reduced controllers of order 18 are obtained. One is by directly applying balanced truncation (BT) (see, e.g. Obinata and Anderson (2001, pp. 7–8)) to the full controller K , and the other is by frequency-weighted balanced truncation (FWBT) (Obinata & Anderson, 2001, pp. 106–107) with the input and output weights being (I − P0 (2, 2)K )−1 P0 (2, 1) and (I − P0 (2, 2)K )−1 P0 (1, 2) respectively (to minimize first order closed loop difference). Fig. 2 shows the magnitude Bode diagram of the error of closed loop systems due to different reduced controllers, as well as the magnitude Bode diagram of the full closed loop system for reference. The left and right figures correspond to the load sensitivity function (from input 1 to output 1) and sensitivity function (from input 1 to output 2) respectively (see Garpinger (2009, pp. 21–22) and Åström and Murray (2009, pp. 316–317) for detail). It is observed that the approximation qualities of the three controller reduction methods are similar in this example. Finally, the integrals of the absolute value of output 1 (i.e. the objective the Youla controller seeks to minimize) due to the full controller, proposed reduced controller and the balanced truncated controllers are 0.9687, 0.9686, 0.9753 and 0.9702 respectively. The lower objective value of the proposed controller is probably due to violation of some design constraint or the inexact (FIR based (Wernrud, 2008)) optimization of the Youla controller. 6.2. Reduction instances from public sources In this subsection some previously studied instances of the controller reduction problem in the setup of Fig. 1 are considered. The data of the plants P are available, in continuous-time state space matrices, from known sources. There are four instances considered, including HIMAT (an aircraft model, data available from the MathWorks website), ROC1, ROC2 and ROC3 from Complib (Leibfritz, 2004). In the experiment, three controller reduction schemes are evaluated. They include the proposed SVD based method, HIFOO (Gumussoy, Henrion, Millstone, & Overton, 2009) and the MATLAB’s balanced truncation methods (taking the best result of ncfmr and balred). ncfmr is an implementation of normalized coprime factors balanced truncation (NCFBT). On the other hand, balred is the basic balanced truncation (BT). Other (weighted) balanced truncation schemes have been reported (Zhou et al., 1995), but only ncfmr and balred are readily available in MATLAB. For the HIFOO test, HIFOO directly designs reduced order controllers from the plants P, without using the information of any full controller. HIFOO is applied five times (with the previous output being the next initial guess) for the HIMAT, ROC1 and ROC2
1074
K.C. Sou, A. Rantzer / Automatica 48 (2012) 1069–1076
(a) Load sensitivity function error magnitude.
(b) Sensitivity function error magnitude.
Fig. 2. Errors of the sensitivity functions in the example in Section 6.1. The sensitivity functions due to the full controller are also shown for reference. Table 1 Closed loop H∞ norm for different controller reduction schemes in the ROC2 example. The full order controller has order 9, with a closed loop H∞ norm of 0.04104. Method
SVD
HIFOO
NCFBT/BT
Order = 4 Order = 6 Order = 7
0.04209 0.04121 0.04104
0.04574 0.04386 0.04240
0.06141 0.04396 0.04105
Table 2 Computational time for different controller reduction schemes in the ROC2 example.
Table 4 Computational time for different controller reduction schemes in the ROC3 example. Method
SVD
HIFOO (s)
NCFBT/BT (s)
Order = 1 Order = 6 Order = 7
N/A 166 s 110 s
702 2921 3383
≈0 ≈0 ≈0
Table 5 Closed loop H∞ norm for different controller reduction schemes in the HIMAT example. The full order controller has order 10, with a closed loop H∞ norm of 0.7885.
Method
SVD (s)
HIFOO (s)
NCFBT/BT (s)
Method
SVD
HIFOO
NCFBT
Order = 4 Order = 6 Order = 7
112 151 144
2043 2964 2800
≈0 ≈0 ≈0
Order = 3 Order = 8
N/A 0.8242
1.045 0.8487
Unstable Unstable
Table 3 Closed loop H∞ norm for different controller reduction schemes in the ROC3 example. The full order controller has order 8, with a closed loop H∞ norm of 45.70. Method
SVD
HIFOO
NCFBT/BT
Order = 1 Order = 6 Order = 7
N/A 52.88 45.71
Unstable 8744 1920
Unstable Unstable 46.78
instances, and ten times for the ROC3 instance. Iterative application of HIFOO is suggested (e.g. Henrion (2006)) because HIFOO employs randomized local searches. For the NCFBT/BT test, full controllers are designed as optimal H∞ norm controllers. Then NCFBT/BT is used to construct reduced controllers. To test the proposed scheme, the plants P are first discretized with a sampling time of 0.01 s (using c2d with the tustin option in MATLAB), and then full order discrete-time optimal controllers are designed. Then the proposed scheme generates discrete-time reduced controllers, which are converted back to continuous-time to evaluate the closed loop H∞ norm. The numerical experiments exhibit advantages, as well as limitations, of the proposed controller method when compared with HIFOO and NCFBT/BT. The ROC2 and ROC3 test cases show the potential of the proposed method. Tables 1 and 2 summarize the closed loop H∞ norm performance and CPU time for the ROC2 case. Tables 3 and 4 analogously summarize the ROC3 case. When compared with HIFOO in these examples, the proposed scheme produces more accurate reduced controllers in a much shorter time (the contrast in the ROC3 case is obvious). When compared with NCFBT/BT, the proposed method generates more
Table 6 Closed loop H∞ norm for different controller reduction schemes in the ROC1 example. The full order controller has order 8, with a closed loop H∞ norm of 1.127. Method Order Order Order Order
=1 =2 =6 =7
SVD
HIFOO
NCFBT
N/A N/A 1.281 1.129
1.426 1.244 1.113 1.196
4.83 35.91 Unstable Unstable
accurate reduced controllers, but balanced truncation requires less time to compute. However, the proposed method fails to find any first order reduced controller for the ROC3 case. This is not surprising though, because the proposed method is a local search in nature. The magnitude Bode diagrams of the closed loop system in the ROC2 example (with reduced controllers of order 4) are shown in Fig. 3. These Bode diagrams depict some sensitivity functions of an aircraft model described in Gangsaas, Bruce, Blight, and Ly (1986). It can be seen that the sensitivity functions due to the proposed reduced controller are the best low order substitutes for the full order ones in this example. On the other hand, the HIMAT and ROC1 cases (in Tables 5 and 6) display the limitation of the proposed method. For example, in the HIMAT case with the current implementation, the proposed method cannot generate reduced controllers of order lower than 8 (the full order being 10). On the contrary, HIFOO manages to find a third order reduced controller with closed loop H∞ norm being 1.0451. For certain applications this may be satisfactory. In addition, in the ROC1 case HIFOO finds acceptable reduced controllers of order one and two, while the proposed
1075
Magnitude
Magnitude
K.C. Sou, A. Rantzer / Automatica 48 (2012) 1069–1076
Frequency (rad/s)
Frequency (rad/s)
(b) Sensitivity function magnitude from vertical gust noise input wn to normal vertical acceleration nz .
Magnitude
Magnitude
(a) Sensitivity function magnitude from horizontal gust input un to normal vertical acceleration nz .
Frequency (rad/s)
Frequency (rad/s)
(c) Sensitivity function magnitude from servo noise input δn to normal vertical acceleration nz .
(d) Sensitivity function magnitude from pitch rate sensor noise qn to normal vertical acceleration nz .
Fig. 3. Sensitivity functions due to the full controller and reduced controllers of order four in the ROC2 example.
method cannot find any. On the other hand, balanced truncation deliver quite lackluster performance in these test cases. In summary, the examples show that the proposed controller reduction scheme manages to find the middle ground between a global optimization based approach (i.e. HIFOO) and local approximation approach (i.e. NCFBT/BT). Hence, it can be useful in complementing the existing schemes. 6.3. Further discussions on algorithmic details Finally, note that the proposed method is applied in Section 6.2 iteratively. First, with the given full controller matrix L, the ‘‘trace heuristic’’ version of problem (18) in Section 5.1 is solved for an initial guess of X . Then, the local update problem in (21) in Section 5.2 is solved several times (e.g. three in this example) for further reduction of the sve(F (X )) (pertaining to the order of the reduced controller). After X is determined, the SVD based ˆ The above process, computation in Section 4 is applied to find L. ˆ from L to L, can be repeated until no more significant progress can be made. The computation time for the SVD method in Tables 2 and 4 include all of the operations above, for all the iterations. It is instructive to see the effect of the local update (21) in Section 5.2, by comparing the achievable reduced controller order with and without the update procedure. Without using the update procedure, for the ROC2 example the reduced controllers generated by the proposed method would have orders of 9, 9, and 8 (as opposed to 4, 6 and 7 in Table 1). For the ROC3 example, the corresponding reduced orders would be 7 and 8 (as opposed to 6 and 7 in Table 3). This demonstrates that computation effort should be spent on solving (18) to as close to optimality as possible.
7. Conclusion This paper demonstrates that it is possible to obtain a closed loop stable controller reduction scheme which linearizes the bounded-real lemma using full controller information—in the spirit of balanced truncation. The proposed controller reduction problem is a minimum rank generalized matrix approximation problem, which can be solved efficiently using SVD. In this paper, numerical examples have demonstrated that proposed controller reduction scheme is on par with the more well-known alternatives such as balanced truncation and HIFOO. Finally, a similar controller reduction setup in H2 norm is also possible. See Sou and Rantzer (2010b) for detail. Acknowledgments Support from the Swedish Research Council through the Linnaeus Centers ACCESS and LCCC are gratefully acknowledged. The authors would also like to thank the anonymous reviewers for helpful comments. Appendix A.1. Proof of Proposition 2 First prove the inequality order(Kˆ )
≤ rank Aˆ k Bˆ k . Bˆ k . The inequality
ˆk , rank A ˆ order(K ) ≤ r holds trivially when r ≥ nˆ k . Therefore, for the rest For convenience, denote r
1076
K.C. Sou, A. Rantzer / Automatica 48 (2012) 1069–1076
of the proof it is assumed that r < nˆ k . Because r < nˆ k , there exists (e.g. by SVD) a unitary matrix S ∈ Rnˆ k ׈nk such that S T S = Inˆ k and S T Aˆ k
Bˆ k = S T Aˆ k
S T Bˆ k =
A˜ k1 0A
B˜ k1 0B
where the zero matrices 0A and 0B have at least nˆ k − r rows. Thus, Kˆ defined by Lˆ is a similarity transform, with transformation matrix S, of a LTI system K˜ , whose system matrices are
A˜ k1 S 0A
,
B˜ k1 0B
, Cˆ k S
ˆ k respectively. Since K˜ has at least nˆ k − r uncontrollable and D states, order(K˜ ) ≤ r. Finally, since Kˆ and K˜ have the same order as they are similarity transforms, it is concluded that order(Kˆ ) ≤ r = rank Aˆ k Bˆ k . The proof of the inequality order(Kˆ ) ≤ rank
ˆ Ak Cˆ k
is omitted because it is similar to the
previous case. Finally, combining the two inequalities leads to the desired result. A.2. Proof of Proposition 3 To show that sve(F (X˜ )) ≤ sve(F (X0 )), denote, for any symmetric matrix M , nne(M ) as the number of nonnegative eigenvalues of M. Then (17) and the discussion immediately before it implies that sve(F (X˜ )) is equal to
A¯ nne ¯ C1
B¯ 1 D11
T
X˜
0
0 ¯ A 1 I C¯ 1
γ
˜ B¯ 1 − X 0 D11
0 γI
(22)
(22) is less than or equal to
nne ci⋆ vi⋆ vi⋆ T +
ci vi vi T +
̸ ⋆ = i ∈I+
j ∈I−
cj vj vj T .
(23)
This relation can be inferred by applying the first constraint in (21) to (17). (23) is less than or equal to
nne λi⋆ vi⋆ vi⋆ T +
̸ ⋆ = i ∈I+
λ i vi vi T +
j ∈I−
λj vj vj T
(24)
and the inequality is true because of the following. X0 , together with the λi ’s (i.e. the eigenvalues of R0 in (20)), is feasible for (21). ̸=⋆ Therefore, for all j ∈ I− , both cj and λj are negative. For all i ∈ I+ , both ci and λi are nonnegative. However, it might be possible that ci⋆ < 0 ≤ λi⋆ , and hence this inequality can be strict. Finally, (24) is equal to sve(F (X0 )) because of (19), (20) and (22). This concludes that sve(F (X˜ )) ≤ sve(F (X0 )). References Ajala, O., Schuler, S., & Allgower, F. (2010). Linf-Gain controller order reduction for discrete-time systems. In American control conference. Åström, K., & Murray, R. (2009). Feedback systems: an introduction for scientists and engineers. Princeton University Press. Anderson, B. D. O., & Liu, Y. (1989). Controller reduction: concepts and approaches. IEEE Transactions on Automatic Control, 34(8). Ankelhed, D., Helmersson, A., & Hansson, A. (2006). Suboptimal model reduction using LMIs with convex constraints, Technical report. Department of Electrical Engineering. Linköping University. Boyd, S., Barratt, C., & Norman, S. (1990). Linear controller design: limits of performance via convex optimization. Proceedings of the IEEE, 78(3), 529–574. Boyd, S., El Ghaoui, L., Feron, E., & Balakrishnan, V. (1994). Linear matrix inequalities in system and control theory. SIAM. Doyle, J., Glover, K., Khargonekar, P., & Francis, B. (1989). State-space solutions to standard H2 and H∞ control problems. IEEE Transactions on Automatic control, 34(8), 831–847.
Enns, D. (1984). Model reduction with balanced realizations: an error bound and a frequency weighted generalization. Conference on Decision and Control. Fazel, M., Hindi, H., & Boyd, S. (2004). Rank minimization and applications in system theory. In American control conference. Boston. Gahinet, P., & Apkarian, P. (1994). A linear matrix inequality approach to H∞ control. International Journal of Robust and Nonlinear Control, 4, 421–448. Gangsaas, D., Bruce, K., Blight, J., & Ly, Uy-Loi (1986). Application of modem synthesis to aircraft control: three case studies. Automatic Control, IEEE Transactions on, 31(11), 995–1014. Garpinger, O. (2009). Design of robust PID controllers with constrained control signal activity. Licentiate Thesis. Reglerteknik. Lund University. Sweden. Goddard, P.J., & Glover, K. (1993). Controller reduction: weights for stability and performance preservation. In Proccedings of 32nd IEEE conference on decision and control (pp. 2903–2908). Gumussoy, S., Henrion, D., Millstone, M., & Overton, M.L. (2009). Multiobjective robust control with HIFOO 2.0. In Proceedings of the IFAC symposium on robust control design. Haifa, Israel. Henrion, D. (2006). Some control design experiments with HIFOO. Technical report 06570. LASS-CNRS. Horn, R., & Johnson, C. (1990). Matrix analysis. Cambridge University Press. Iwasaki, T., & Skelton, R. E. (1994). All controllers for the general H∞ control problem: LMI existence conditions and state space formulas. Automatica, 30(8), 1307–1317. Iwasaki, T., & Skelton, R. E. (1995a). The XY -centring algorithm for the dual LMI problem: a new approach to fixed-order control design. International Journal of Control, 62(6), 1257–1272. Iwasaki, T., & Skelton, R. E. (1995b). A unified approach to fixed-order controller design via linear matrix inequalities. Mathematical Problems in Engineering, 1(1), 59–75. Leibfritz, F. (2004). COMPleib: COnstraint matrix-optimization problem library—a collection of test examples for nonlinear semidefinite programs. Control system design and related problems. Technical report. University of Trier. Department of Mathematics. data available from http://www.complib.de/. Obinata, G., & Anderson, B. D. O. (2001). Model reduction for control system design. Springer. Sandberg, H. (2010). An extension to balanced truncation with application to structured model reduction. IEEE Transactions on Automatic Control, 55(4), 1038–1043. Sou, K.C., & Rantzer, A. (2010a). A singular value decomposition based closed loop stability preserving controller reduction method. American Control Conference, Baltimore, 2010, pp. 1079–1084. Sou, K.C., & Rantzer, A. (2010b). Controller reduction with closed loop error guarantee. 49th Conference on Decision and Control, Atlanta, 2010. Sou, K.C., & Rantzer, A. (2010c). On the minimum rank of a generalized matrix approximation problem in the maximum singular value norm. In 19th international symposium on mathematical theory of networks and systems. Sturm, J. F. (1999). Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones. Optimization Methods and Software, 11–12, 625–653. Wernrud, A. (2008). QTool 0.1 — Reference manual. Department of automatic control. Lund University, Sweden. Löfberg, J. (2004). Yalmip: a toolbox for modeling and optimization in MATLAB. In Proceedings of the CACSD Conference. Taipei, Taiwan. Zhou, K. (1993). Frequency weighted model reduction with L∞ error bounds. Systems and Control Letters, 21, 115–125. Zhou, K., Doyle, J., & Glover, K. (1995). Robust and optimal control. Prentice Hall.
Kin Cheong Sou received a Ph.D. degree in Electrical Engineering and Computer Science at Massachusetts Institute of Technology in 2008. From 2008 to 2010 he was a postdoctoral researcher at Lund University, Lund, Sweden, and since September 2010 he has been a postdoctoral researcher at KTH Royal Institute of Technology, Stockholm, Sweden. His research interests include power system cybersecurity analysis, environment aware building and community, convex/non-convex optimization and model reduction for dynamical systems. Anders Rantzer received a Ph.D. degree in Optimization and Systems Theory in 1991 from the Royal Institute of Technology (KTH), Stockholm, Sweden. After postdoctoral positions at KTH and at IMA, University of Minnesota, he joined the Department of Automatic Control at Lund University in 1993. He was appointed professor of Automatic Control in Lund 1999. The academic year of 2004/05 he was a visiting associate faculty member at Caltech. Since 2008 he has coordinated the Linnaeus center LCCC at Lund University. Rantzer has been serving as an associate editor of IEEE Transactions on Automatic Control and several other journals. He is a winner of the SIAM Student Paper Competition, the IFAC Congress Young Author Price and the IET Premium Award for the best article in IEE Proceedings — Control Theory & Applications during 2006. He is a Fellow of IEEE and a member of the Royal Swedish Academy of Engineering Sciences. His research interests are in modeling, analysis and synthesis of control systems, with particular attention to uncertainty, optimization and distributed control.