- Email: [email protected]
Generalized pole placement via static output feedback: A methodology based on projections
Automatica 42 (2006) 2143 – 2150 www.elsevier.com/locate/automatica
Brief paper
Generalized pole placement via static output feedback: A methodology based on projections夡 Kaiyang Yang a,b,∗ , Robert Orsi a a Research School of Information Sciences and Engineering, The Australian National University, Canberra ACT 0200, Australia b National ICT Australia Limited, Locked Bag 8001, Canberra ACT 2601, Australia
Received 30 September 2005; received in revised form 14 March 2006; accepted 29 June 2006 Available online 30 August 2006
Abstract This paper presents an algorithm for solving static output feedback pole placement problems of the following rather general form: given n subsets of the complex plane, find a static output feedback that places in each of these subsets a pole of the closed-loop system. The algorithm presented is iterative in nature and is based on alternating projection ideas. Each iteration of the algorithm involves a Schur matrix decomposition, a standard least-squares problem and a combinatorial least-squares problem. While the algorithm is not guaranteed to always find a solution, computational results are presented demonstrating the effectiveness of the algorithm. 䉷 2006 Elsevier Ltd. All rights reserved. Keywords: Static output feedback; Pole placement; Stabilization; Alternating projections
1. Introduction
In this paper we will actually consider the following generalized static output feedback pole placement problem.
There has been a great deal of research done on the problems of pole placement and stabilization via static output feedback. An overview of theoretical results, existing algorithms and historical developments can be found in Byrnes (1989), Rosenthal and Willems (1998), Syrmos, Abdallah, Dorato, and Grigoriadis (1997), de Oliveira and Geromel (1997), Rosenthal and Sottile (1998) and Eremenko and Gabrielov (2002). Our main interest here is in algorithms, and in this regard, for pole placement, the survey paper Rosenthal and Willems (1998) states that existing sufficiency conditions are mainly theoretical in nature and that there are no good numerical algorithms available in many cases when a problem is known to be solvable. Despite the great deal of work that has been done in this area, new algorithms for these important problems are still of interest.
Problem 1. Given A ∈ Rn×n , B ∈ Rn×m , C ∈ Rp×n and closed subsets C1 , . . . , Cn ⊂ C, find K ∈ Rm×p such that i (A + BKC) ∈ Ci for i = 1, . . . , n.
夡 A preliminary version of this paper was presented at the 16th IFAC World Congress. This paper was recommended for publication in revised form by Associate Editor Yasumasa Fujisaki under the direction of Editor Roberto Tempo. ∗ Corresponding author. Research School of Information Sciences and Engineering, The Australian National University, Canberra ACT 0200, Australia. E-mail addresses: [email protected] (K. Yang), [email protected] (R. Orsi).
0005-1098/$ - see front matter 䉷 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.automatica.2006.06.021
Here i (A + BKC) denotes the ith eigenvalue of A + BKC. Problem 1 encompasses many types of pole placement problems. Here are some examples: (1) Classical pole placement: Ci = {ci },
ci ∈ C.
(2) Stabilization type problems for continuous time and discrete time systems: C1 = · · · = Cn = {z ∈ C | Re(z) − },
>0
and C1 = · · · = Cn = {z ∈ C | |z| }, respectively. (3) Relaxed classical pole placement: Ci = {z ∈ C | |z − ci | ri }.
0 < < 1,
(1)
2144
K. Yang, R. Orsi / Automatica 42 (2006) 2143 – 2150
Fig. 1. A hybrid problem.
Here each region Ci is a disk centered at ci ∈ C with radius ri 0. (4) Hybrid problems: For example, problems of the type shown in Fig. 1. Here c ∈ C and the aim is to place a pair of poles at c and c, ¯ and to place the remaining poles in the truncated cone C: C1 = {c},
C2 = {c}, ¯
C3 = · · · = Cn = C.
As far as the authors are aware, pole placement in the generality presented in Problem 1 has not previously been considered. The most closely related results from the existing literature are as follows. Early work in Gutman and Jury (1981) considers pole placement in a single region specified by polynomials. While a Lyapunov-type necessary and sufficient condition is given for a matrix to have its eigenvalues in such a region, this condition is polynomial in the matrix in question and hence not readily amenable to controller design. In Chilali and Gahinet (1996), LMI conditions are presented that are sufficient though not necessary for pole placement in various convex regions. These results cover state feedback and fullorder dynamic output feedback but not static output feedback. Another LMI-based approach, again for a single, though this time, possibly disconnected region, is considered in Bosche, Bachelier, and Mehdi (2004). A method for placing poles in distinct convex regions (each region is specified using linear programming constraints or second-order cone constraints) is given in Hassibi, How, and Boyd (1999); however, the method is based on eigenvalue perturbation results and hence appears largely limited to cases where the open-loop poles are already quite close to the desired poles. In Satoh, Okubo, and Sugimoto (2003), pole placement in distinct convex regions (each region is a disk or a half plane) is achieved via a rank constrained LMI approach though the results are only for state feedback. This paper presents an algorithm for Problem 1. The approach employed here is quite different to each of the approaches mentioned above. Problem 1 is shown to be equivalent to finding a point in the intersection of two particular sets, one of which is a simple convex set and the other a rather complicated nonconvex set. The algorithm is iterative in nature and is based on an alternating projection like scheme between these two sets. Each iteration of the algorithm involves a Schur matrix decomposition and a standard least-squares problem. If the Ci ’s are not all equal, each iteration also requires a combinatorial least-squares matching step.
Alternating projection type ideas have been employed previously for output feedback stabilization, see in particular Grigoriadis and Skelton (1996), Grigoriadis and Beran (1999) and Orsi, Helmke, and Moore (2006).1 A distinguishing feature of our algorithm is that, unlike these methods, our algorithm does not involve LMIs. (A further technical difference is that rather than solving feasibility problems that involve symmetric matrices, the problem solved by the algorithm is a feasibility problem involving nonsymmetric matrices.) The algorithm can be applied to problems with rather general choices for the Ci regions. In fact the only formal requirement is the following, which we state in the form of an assumption. Assumption 1. It is possible to calculate projections onto each of the Ci ’s: given z ∈ C it is possible to find zi ∈ Ci such that |z − zi | |z − c| for all c ∈ Ci . In particular, the Ci ’s must be closed sets though they need not be convex or even connected. For given A, B, C and Ci ’s, Problem 1 may or may not have a solution. Indeed, one would expect that determining whether a particular instance of Problem 1 is solvable is in general difficult. For example, the problem of determining whether the classical pole placement problem is solvable for particular A, B, C and desired poles has recently been shown to be NPhard (Fu, 2004). Given the difficulty of Problem 1, an efficient (i.e., polynomial time) algorithm that is able to correctly solve all instances of the problem cannot be expected, and while the algorithm presented here is often quite effective in practice, it is not guaranteed to find a solution even if a solution exists. The structure of the paper is as follows. The remainder of this section lists some notation which is used in the rest of the paper. Projections play a key part in the algorithm and Section 2 contains general properties of projections and recalls how alternating projections can be used to find a point in the intersection of a finite number of closed (convex) sets. Section 3 presents the solution methodology. To motivate the solution methodology we first restrict our attention to systems that have a symmetric state space representation (A = AT , C = B T ) and present an algorithm for this easier class of problems. The general problem is then considered. Section 4 contains computational results of applying the algorithm to various instances of Problem 1. The paper ends with some concluding remarks. Notation: Sr is the set of real symmetric r × r matrices. diag(v) for v ∈ Cr denotes the r × r diagonal matrix whose ith diagonal term is vi . For Z ∈ Cr×s , vec(Z) ∈ Crs consists of the columns of Z stacked below each other. Y ⊗ Z denotes the Kronecker product of Y and Z. For Z ∈ Cr×s , Re(Z) ∈ Rr×s and Im(Z) ∈ Rr×s denote, respectively, the real and imaginary parts of Z. 1 They were first used in control design to solve certain convex problems, see Grigoriadis and Skelton (1994).
K. Yang, R. Orsi / Automatica 42 (2006) 2143 – 2150
2. Projections Projections play a key part in the algorithm. This section contains some general properties of projections that are used throughout the paper. Let x be an element in a Hilbert space H and let D be a closed (possibly nonconvex) subset of H. Any d0 ∈ D such that x − d0 x − d for all d ∈ D will be called a projection of x onto D. In the cases of interest here, namely that H is a finite dimensional Hilbert space, there is always at least one such point for each x. If D is convex as well as closed then each x has exactly one such minimum distance point. A function PD : H → H will be called a projection operator (for D) if for each x ∈ H , PD (x) is a projection of x onto D. In this paper the problems we are interested in solving are feasibility problems of the following general form. Problem 2. Given closed sets D1 , . . . , DN in a finitedimensional Hilbert space H, find a point in the intersection N i=1 Di (assuming the intersection is non-empty). (In fact we will only be interested in the case N = 2.) When all the Di ’s in Problem 2 are convex, a point in the intersection of the Di ’s can be found via alternating projections. Theorem 3. Let D1 , . . . , DN be closed convex sets in a finite dimensional Hilbert space H. Given constants t1 , . . . , tN in the interval (0, 2), for i = 1, . . . , N, define the operators Ri = (1 − ti )I d + ti PDi . If N i=1 Di is non-empty, then starting from an arbitrary initial value, the sequence defined by xi+1 = Ri (mod N)+1 (xi ) converges to an element in N D . i i=1
projections onto these sets are computable, a solution method is to alternately project onto D1 and D2 . Theorem 4 ensures that the distance xn − yn is nonincreasing with n. While the convergence of xn −yn to 0 is not guaranteed, if this quantity does converge to 0 then a solution to the problem will have been found. 3. Methodology 3.1. The symmetric problem Systems with a symmetric state space realization, that is, systems with state space matrices satisfying A = AT , C = B T , occur in various contexts, for example RC-networks. In order to motivate our solution methodology, we first consider the following special case of Problem 1 for symmetric systems.2 Problem 5. Given A ∈ Sn , B ∈ Rn×m , and closed subsets C1 , . . . , Cn ⊂ R, find K ∈ Sm such that i (A + BKB T ) ∈ Ci for i = 1, . . . , n. Note that in Problem 5 the static output feedback matrix K is required to be symmetric and that the Ci ’s are assumed to be subsets of the real numbers. This latter assumption is natural as given any symmetric K, A + BKB T is also symmetric and hence has real eigenvalues. We assume A, B and C1 , . . . , Cn are given and fixed. Let L = {Z ∈ Sn | Z = A + BKB T for some K ∈ Sm } and let M denote the set of symmetric matrices with eigenvalues in the specified regions C1 , . . . , Cn ,
Proof. See Gubin, Polyak, and Raik (1967). An alternate proof can also be found in Youla and Webb (1982).
M = {Z ∈ Sn | i (Z) ∈ Ci , i = 1, . . . , n}.
Here Id denotes the identity operator on H; I d(x) = x for all x ∈ H . Each Ri of the theorem is a projection operator onto Di if and only if ti = 1. If ti = 1 then Ri (x) is referred to as a relaxed projection of x. Theorem 3 is a generalization of the classical method of alternating projections which corresponds to the case t1 = · · · = tN = 1. We will see later how the freedom to choose the ti ’s not equal to 1 can be useful. When one or more Di ’s are nonconvex, Theorem 3 no longer applies and starting the algorithm of Theorem 3 from certain initial values may result in a sequence of points that does not converge to a solution of the problem (Combettes & Trussell, 1990). If, however, there are only two sets, then the following distance reduction property always holds (Orsi & Yang, 2006).
The symmetric problem can be stated as follows:
Theorem 4. Let D1 and D2 be closed, nonempty, possibly nonconvex sets in a finite dimensional Hilbert space H, and let PD1 and PD2 be associated projection operators. For any initial value y0 ∈ H and for all n 0, let xn+1 = PD1 (yn ) and yn+1 =PD2 (xn+1 ). Then for all n 1, xn+1 −yn+1 xn+1 − yn xn − yn . Suppose one is interested in solving Problem 2 in the case of two sets, D1 and D2 , when one or both sets are nonconvex. If
2145
Find X ∈ L ∩ M. We now show that projections onto both sets L and M can be calculated and hence that alternating projections can be employed as a solution method. From now on, Sn will be regarded as a Hilbert space with inner product Y, Z = tr(Y Z) and associated norm Z = Z, Z 1/2 (the Frobenius norm). For i = 1, . . . , n, PCi will denote a projection operator for Ci . The set L is an affine subspace and hence is convex. Projection of X ∈ Sn onto L involves solving a standard leastsquares problem. The following result is fairly straightforward to establish. Lemma 6. The projection of X ∈ Sn onto L is given by PL (X) = A + BKB T where K ∈ Sm is a minimizer of (B ⊗ B)vec(K) − vec(X − A)2 . ( · 2 denotes the standard vector 2 norm.) 2 See Mahony and Helmke (1995) for more regarding the classical pole placement problem for symmetric systems and a rather surprising result regarding arbitrary pole placement for such systems.
2146
K. Yang, R. Orsi / Automatica 42 (2006) 2143 – 2150
While the set M is not convex, it is still possible to calculate projections onto M. How to do this is shown in Theorem 8. The result is based on the following result of Hoffman and Wielandt. Lemma 7. Suppose Y, Z ∈ Sn have eigenvalue-eigenvector decompositions Y = V DV T , D = diag(Y1 , . . . , Yn ), and Z = Z n×n W EW T , E = diag(Z are or1 , . . . , n ), where V , W ∈ R Y Y Z Z thogonal and 1 · · · n and 1 · · · n . Then D − EY − Z,
(2)
where · denotes the Frobenius norm.
Theorem 8. Given Y ∈ Sn , let Y = V DV T be an eigenvalueeigenvector decomposition of Y with D = diag(1 , . . . , n ). Let be a permutation of {1, . . . , n}such that, amongst all possible permutations, it minimizes nk=1 |k − PC(k) (k )|2 . ˆ T where Dˆ = diag(PC (1 ), . . . , Define PM (V , D) = V DV (1) PC(n) (n )). Then PM (V , D) is a best approximant in M to Y in the Frobenius norm. Proof. Let Y be as in the theorem statement. As PM (V , D) ∈ M, it remains to show for all Z ∈ M.
(3)
Without loss of generality, suppose the eigenvalues of Y are ordered, i.e., 1 · · · n . Similarly, for Z ∈ M, let Z = W EW T be an eigenvalue-eigenvector decomposition with E = Z Z Z diag(Z 1 , . . . , n ) and 1 · · · n . As the Frobenius norm is orthogonally invariant, ˆ Y − PM (V , D) = D − D.
(4)
Let be a permutation of {1, . . . , n} such that Z k ∈ C(k) , k = 1, . . . , n. (Such a permutation exists as Z ∈ M and hence must have an eigenvalue in each of the Ci ’s.) Define Eˆ = diag(PC(1) (1 ), . . . , PC(n) (n )). It follows from the definition of Dˆ that ˆ D − E. ˆ D − D
(5)
As for each k, |k − PC(k) (k )||k − Z k |, it also follows that ˆ D − ED − E.
L = {Z ∈ Rn×n | Z = A + BKC for some K ∈ Rm×p }, and re-define M to be the set of complex matrices with eigenvalues in the specified regions C1 , . . . , Cn , M = {Z ∈ Cn×n | i (Z) ∈ Ci , i = 1, . . . , n}. Problem 1 can now be stated as
Proof. See for example Horn and Johnson (1985, Corollary 6.3.8).
Y − PM (V , D) Y − Z
From now on Cn×n will be regarded as a Hilbert space with inner product Y, Z = tr(Y Z ∗ ). In this subsection we re-define L to be the set of all possible closed-loop system matrices,
(6)
Combining (4)–(6) and inequality (2) from Lemma 7 gives the inequality in (3) and the proof is complete. Note that to calculate PM (V , D) we keep the original orthogonal matrix V and simply modify the diagonal matrix D to ˆ The fact that V remains unchanged motivates our solution D. method for the general nonsymmetric case. 3.2. The general nonsymmetric problem Consider again Problem 1. Throughout this subsection it is assumed A, B, C and C1 , . . . , Cn are given and fixed.
Find X ∈ L ∩ M. A solution strategy to solve Problem 1 would be to employ an alternating projection scheme, alternatively projecting between L and M. A difficulty occurs in trying to do this. While L is convex and M is nonconvex, just as for the symmetric problem, M in this case is a substantially more complicated set and how to calculate projections onto this set is a difficult unsolved problem. As will be verified by the experiments, an alternating projection like scheme can still be quite successful if for M a suitable substitute mapping is used in place of an actual projection operator. Before proceeding, recall Schur’s result (Horn & Johnson, 1985). Theorem 9. Given Z ∈ Cn×n with eigenvalues 1 , . . . , n in any prescribed order, there is a unitary matrix V ∈ Cn×n and an upper triangular matrix T ∈ Cn×n such that Z = V T V ∗ , and Tkk = k for k = 1, . . . , n. The following map is proposed as a substitute for a projection map onto M. Though it is not a true projection map, the notation PM will still be used. Our choice for PM is motivated by the projection operator for M for the symmetric problem and the idea for using such a substitute mapping originates in Orsi and Yang (2006) (see also Orsi (2006)). We stress that as PM is not a true projection map, unlike for the symmetric problem, the distance reduction property of Theorem 4 may not hold. While PM has various desirable properties (see below), the theoretical convergence properties of the algorithm are currently unclear, providing interesting questions for future research. In what follows, Assumption 1 is assumed to hold. Definition 10. Suppose V ∈ Cn×n is unitary and T ∈ Cn×n is upper triangular. Let be a permutation of {1, . . . , n} such that amongst all possible permutations, it minimizes n
|Tkk − PC(k) (Tkk )|2 .
(7)
k=1
Define PM (V , T ) = V Tˆ V ∗ where Tˆ is upper triangular and given by PC(k) (Tkk ) if k = l, Tˆkl = Tkl otherwise.
K. Yang, R. Orsi / Automatica 42 (2006) 2143 – 2150
Note that PM maps into the set M. Note also that, just as for the symmetric problem, finding involves solving a combinatorial least-squares problem. (This will be discussed further later in the section.) In order to apply PM to Z ∈ Cn×n , a Schur decomposition of Z must first be found. A given Z may have a nonunique Schur decomposition and Z = V1 T1 V1∗ = V2 T2 V2∗ does not necessarily imply PM (V1 , T1 ) = PM (V2 , T2 ). Hence, PM may give different points for different Schur decompositions of the same matrix. (Though it was not discussed at the time, similar comments apply to the PM mapping for the symmetric case.) This is not so important as different Schur decompositions lead to points in M of equal distance to the original matrix. Also note that of all the points in M that have a Schur decomposition with the same V matrix, PM (V , T ) is closest (or at least equal closest) to the original point Z = V T V ∗ . For a proof of these last two claims, see Orsi and Yang (2006) or Orsi (2006). As for the symmetric problem, projection of X ∈ Cn×n onto L involves solving a standard least-squares problem. Lemma 11. The projection of X ∈ Cn×n onto L is given by PL (X) = A + BKC where K ∈ Rm×p is a minimizer of (C T ⊗ B)vec(K) − vec(Re(X) − A)2 . Here is our algorithm for Problem 1. Algorithm. Problem Data: A ∈ Rn×n , B ∈ Rn×m , C ∈ Rp×n and C1 , . . . , Cn ⊂ C. Initialization: Choose a randomly generated matrix Y ∈ Rn×n . For example, draw each entry of Y from a normal distribution of zero mean and variance 1. repeat 1. X := PL (Y ). 2. Calculate a Schur decomposition of X: X = V T V ∗ . 3. Y := PM (V , T ). until X − Y < . see Definition 10, Note that as Y = PM (V , T ) = V Tˆ V ∗ , X −Y =V T V ∗ −V Tˆ V ∗ =T − Tˆ =( |Tkk − Tˆkk |2 )1/2 . As the Tkk ’s are the eigenvalues of X ∈ L and the Tˆkk ’s are the eigenvalues of Y ∈ M, the algorithm stops when X (which equals A+BKC for some K) has eigenvalues sufficiently close to those of a matrix that satisfies the pole placement constraints. In particular, each eigenvalue of such an X cannot violate the pole placement constraints by more than . As mentioned previously, PM involves finding a permutation that minimizes (7). The first step in solving this combinatorial least-squares problem is calculating the squared distance of each Tkk to each subset Cl and placing these values in an n × n cost matrix D: Dkl := |Tkk − PCl (Tkk )|2 .
(8)
The problem is now equivalent to finding a permutation such that Dk (k) is minimal.
2147
The problem of finding a minimizing given a cost matrix D is a linear assignment problem which can be solved in O(n3 ) time using the so-called Hungarian method, see Martello and Toth (1987) for details. Note that if the Ci ’s are not all distinct, for example this occurs for stabilization problems and what here have been called hybrid problems, the complexity of the matching problem is reduced. In fact for stabilization problems, all the Ci ’s are the same and no matching step is required. For the hybrid problem shown in Fig. 1, it is only necessary to check n(n − 1) possibilities corresponding to which two Tkk ’s are matched to c and c. ¯ Hence for this hybrid problem, the direct approach is faster than using the Hungarian method. Also note that, given a cost matrix D, an alternative to the Hungarian method is the following faster suboptimal matching strategy. Find the (or a) smallest entry in D, say Dk˜ l˜. Match Tk˜ k˜ with Cl˜ and cross out row k˜ and column l˜ of D. Now only consider the uncrossed out entries in D and repeat, until all n matches have been made. This method does not always find the optimal matching though it can often be quite an effective substitute for the Hungarian method. It will be termed suboptimal matching. Surprisingly, as will be shown in the next section, by using suboptimal matching in the algorithm it was possible to solve a particular problem which was not solvable using the Hungarian method. 4. Computational results This section contains computational results of applying the algorithm to various instances of Problem 1. We include results for classical pole placement, continuous time stabilization, and a hybrid problem. (Additional results can be found in Yang & Orsi (2005) and Yang, Orsi, & Moore (2004).) The algorithms for each problem were implemented in Matlab 6.5 and all results were obtained using a 3.06 GHz Pentium 4 machine. Throughout this section a randomly generated matrix will be a matrix whose entries are drawn from a normal distribution of zero mean and variance 1. 4.1. Classical pole placement: random problems This subsection contains results for some randomly generated classical pole placement problems. A 1000 problems with n=6, m = 4 and p = 3 were created. Each problem was created as follows. A, B, C and K˜ were generated randomly. A scalar multiple of the identity was added to A to ensure the stability ˜ was equal to = 0.1. The desired poles degree of A + B KC ˜ were taken to be the poles of A + B KC. An attempt was made to solve each problem using up to 10 different initial conditions and a maximum of a 1000 iterations per initial condition. With the termination parameter set to = 10−3 , the overall success rate was 91%. (The success rate increases to 95% if up to 20 initial conditions are used.) For the problems that were solved, the average number of iterations taken was 1.8 × 103 and the average time taken was 1.5 CPU seconds.
2148
K. Yang, R. Orsi / Automatica 42 (2006) 2143 – 2150
Table 1 A comparison of performance for different n, m and p. i denotes the number of iterations and listed are the number of problems that converged within different iteration ranges. ‘NC’ denotes the number of problems that did not converge within 1000 iterations. T denotes the average convergence time in CPU seconds n m p
3 1 1
6 1 1
9 1 1
3 2 1
6 4 3
9 5 4
1 i 10 10 < i 100 100 < i 1000 NC
996 3 1 0
990 7 3 0
991 8 1 0
993 3 4 0
989 11 0 0
982 17 1 0
T
0.0007
0.0007
0.001
0.0006
0.001
0.002
4.2. Classical pole placement: a particular problem from the literature The following problem is taken from Sridhar and Lindorff (1973) and is of interest as the set of desired poles overlaps with the set of open-loop poles. The system matrices are the following: ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ 1 0 0 0 1 0 1 0 T 0 ⎥ ⎢0 2 0 ⎢0 1⎥ ⎢1 0⎥ A=⎣ ⎦, B = ⎣ ⎦, C = ⎣ ⎦ . 0 0 −3 0 1 0 0 1 1 1 0 1 0 0 0 −4 The aim is to place the closed-loop poles at {−1, −2, −3, −5}. While initial attempts to solve this problem failed, the problem was solved by using suboptimal matching and replacing Step 3 of the algorithm with the relaxed projection ‘Y := (1 − t)X + tP M (V , T )’, t ∈ (0, 2) constant. Strictly speaking, this is not a relaxed projection in the true sense as PM (V , T ) may not be a projection of X = V T V ∗ , however, the idea is clearly the same. Solutions were successfully found by taking t close to 0; t =0.1, 0.2 and 0.3 can all be used to successfully find a solution. Likelihood of success increases and speed of convergence decreases the closer t is to 0. With t = 0.3 and = 10−3 , solutions can typically be found in about 1.2 × 104 iterations and about 4.7 CPU seconds. Note: When employing relaxed projections, the loop termination criterion of the algorithm should be replaced by ‘until X − PM (V , T ) < ’ as now X − Y = X − ((1 − t)X + tP M (V , T )) = tX − PM (V , T ). 4.3. Continuous time stabilization: random problems In the comparison paper de Oliveira and Geromel (1997), a number of methods for continuous time stabilization via static output feedback were compared. In this subsection we repeat the same numerical experiments of de Oliveira and Geromel (1997) using our algorithm, to see how it compares. One thousand problems were generated for each of a number of different choices of system dimensions n, m and p. In each case, A, B, C and K˜ were generated randomly. A scalar multiple of the identity was added to A to ensure it had a stability degree ˜ All stable A’s were of one. A was then replaced by A − B KC. discarded.
In applying our algorithm, the Ci ’s were chosen as in (1) with = 0.1. As closed-loop systems were only required to be stable, the algorithm was terminated as soon as all poles had real part less than zero. An attempt was made to solve each problem using up to 10 different initial conditions and a maximum of a 100 iterations per initial condition. Results are given in Table 1. As can be seen, all problems were solved. Most problems were solved in 10 iterations or less and average solution times were very small. These results are as good as the results for the best two algorithms tested in de Oliveira and Geromel (1997). Other values for , such as 0.2, 0.5, 1 and 2, were also tried and produced just as good results, indicating a certain robustness of the algorithm with respect to this parameter. While the above results are very encouraging, the algorithm may not work as well for all problems. For example, in Yang et al. (2004) the method of problem generation is different and the algorithm does not perform as well. 4.4. A hybrid problem So far we have presented results for two traditional classes of problems. In this subsection we demonstrate the generality of the algorithm by considering a less standard problem, namely a hybrid problem of the type shown in Fig. 1. The problem parameters are c = −0.5 + i3, C = {z ∈ C | Re(z) − 2 and |Im(z)| |Re(z)|}, n = 13, m = 3 and p = 5. To ensure solvability, B, C and K˜ were ran˜ with domly generated and A set to A = V˜ T˜ V˜ T − B KC ˜ ˜ V a random orthogonal matrix and T a real block upper triangular matrix with a spectrum satisfying the constraints. (T˜ was assigned the spectrum {−0.5 ± i3, −2, −2 ± i, −2.3, −2.5, −3 ± i3, −3.5 ± i3.1, −4 ± i4} and was created by choosing appropriate 1 × 1 and 2 × 2 blocks for its diagonal. The remaining upper triangular entries of T˜ were chosen randomly.) With = 10−3 , 64% of initial conditions tested resulted in a solution of this problem within 5000 iterations. For the initial conditions that lead to convergence, the average number of iterations taken was 8.6 × 102 and the average time taken was 3.1 CPU seconds. The closed-loop poles corresponding to a particular solution are shown in Fig. 2.
K. Yang, R. Orsi / Automatica 42 (2006) 2143 – 2150 5 4 3 2 1 0 -1 -2 -3 -4 -5
-5
-4.5
-4
-3.5
-3
-2.5
-2
-1.5
-1
-0.5
0
Fig. 2. The closed-loop poles corresponding to a solution for the considered hybrid problem.
Note: The least-squares matching steps were done directly rather than using the Hungarian method. 5. Conclusion In this paper a new methodology for solving a broad class of output feedback pole placement problems has been presented. While the methodology is not guaranteed to find a solution, numerical experiments presented in the paper demonstrate it can be quite effective in practice. A particular strength of the algorithm is that, in addition to being able to solve classical pole placement and stabilization problems, it can also be used to solve less standard pole placement problems. Acknowledgments The authors thank Iven Mareels and Robert Mahony for their helpful guidance to various parts of the pole placement literature. The second author acknowledges the support of the Australian Research Council through Grant DP0450539. National ICT Australia is funded through the Australian Government’s Backing Australia’s Ability initiative, in part through the Australian Research Council. References Bosche, J., Bachelier, O., & Mehdi, D. (2004). Robust pole placement by static output feedback. In Proceedings of the 43rd IEEE conference on decision and control (pp. 869–874). Paradise Island, Bahamas. Byrnes, C. I. (1989). Pole assignment by output feedback. In Nijmeijer, & Schumacher (Eds.), Three decades of mathematical system theory, Lecture notes in control and information science Vol. 135, (pp. 31–78). Berlin: Springer. Chilali, M., & Gahinet, P. (1996). H∞ design with pole placement constraints: An LMI approach. IEEE Transactions on Automatic Control, 41(3), 358–367. Combettes, P. L., & Trussell, H. J. (1990). Method of successive projections for finding a common point of sets in metric spaces. Journal of Optimization Theory and Applications, 67, 487–507.
2149
de Oliveira, M. C., & Geromel, J. C. (1997). Numerical comparison of output feedback design methods. In Proceedings of American control conference (pp. 72–76). New Mexico, USA. Eremenko, A., & Gabrielov, A. (2002). Pole placement by static output feedback for generic linear systems. SIAM Journal of Control Optimization, 41(1), 303–312. Fu, M. (2004). Pole placement via static output feedback is NP-hard. IEEE Transactions on Automatic Control, 49(5), 855–857. Grigoriadis, K. M., & Beran, E. B. (1999). Alternating projection algorithms for linear matrix inequalities problems with rank constraints. In L. El Ghaoui, & S.-I. Niculescu (Eds.), Advances on linear matrix inequality methods in control (pp. 251–267). Philadelphia, PA: SIAM. Grigoriadis, K. M., & Skelton, R. E. (1994). Alternating convex projection methods for covariance control design. International Journal of Control, 60(6), 1083–1106. Grigoriadis, K. M., & Skelton, R. E. (1996). Low-order control design for LMI problems using alternating projections. Automatica, 32(8), 1117–1125. Gubin, L. G., Polyak, B. T., & Raik, E. V. (1967). The method of projections for finding the common point of convex sets. USSR Computational Mathematics and Mathematical Physics, 7, 1–24. Gutman, S., & Jury, E. (1981). A general theory for matrix root-clustering in subregions of the complex plane. IEEE Transactions on Automatic Control, 26(4), 853–863. Hassibi, A., How, J., & Boyd, S. (1999). A path-following method for solving BMI problems in control. In Proceedings of the American control conference (Vol. 2) (pp. 1385–1389), San Diego, USA. Horn, R. A., & Johnson, C. R. (1985). Matrix analysis. Cambridge: Cambridge University Press. Mahony, R. E., & Helmke, U. (1995). System assignment and pole placement for symmetric realisations. Journal of Mathematical Systems, Estimation and Control, 5(2), 1–32. Martello, S., & Toth, P. (1987). Linear assignment problems. Annals of Discrete Mathematics, 31, 259–282. Orsi, R. (2006). Numerical methods for solving inverse eigenvalue problems for nonnegative matrices. SIAM Journal on Matrix Analysis and Applications, 28(1), 190–212. Orsi, R., Helmke, U., & Moore, J. B. (2006). A Newton-like method for solving rank constrained linear matrix inequalities. Automatica, 42 (11). Orsi, R., & Yang, K. (2006). Numerical methods for solving inverse eigenvalue problems for nonnegative matrices. In Proceedings of the 17th international symposium on mathematical theory of networks and systems (pp. 2284–2290). Kyoto, Japan. Rosenthal, J., & Sottile, F. (1998). Some remarks on real and complex output feedback. Systems and Control Letters, 33, 73–80. Rosenthal, J., & Willems, J. C. (1998). Open problems in the area of pole placement. In V.D. Blondel, E.D. Sontag, M. Vidyasager, & J.C. Willems (Eds.), Open problems in mathematical systems and control theory (pp. 181–191). Berlin: Springer. Satoh, A., Okubo, J., & Sugimoto, K. (2003). Transient response shaping in H∞ control by eigenstructure assignment to convex regions. In Proceedings of the 42nd IEEE conference on decision and control (pp. 2288–2293). Maui, USA. Sridhar, B., & Lindorff, D. P. (1973). Pole placement with constraint gain output feedback. International Journal of Control, 18(5), 993–1003. Syrmos, V. L., Abdallah, C. T., Dorato, P., & Grigoriadis, K. (1997). Static output feedack—A survey. Automatica, 33(2), 125–137. Yang, K., & Orsi, R. (2005). Pole placement via output feedback: A methodology based on projections. In Proceedings of the 16th IFAC world congress, Prague, Czech Republic. Yang, K., Orsi, R., & Moore, J. B. (2004). A projective algorithm for static output feedback stabilization. In Proceedings of the second IFAC symposium on system, structure and control (pp. 263–268). Oaxaca, Mexico. Youla, D. C., & Webb, H. (1982). Image restoration by the method of convex projections: Part 1—Theory. IEEE Transactions on Medical Imaging, MI1(2), 81–94.
2150
K. Yang, R. Orsi / Automatica 42 (2006) 2143 – 2150 Kaiyang Yang was born in Shenyang, PR China, in 1979. She received her B.E. degree in Control Science and Engineering from Zhejiang University, PR China, in 2002. From July 2002 to June 2003, she worked for Jinboyuan System Integration as a system engineer. Ms. Yang started her Ph.D. studies at the Australian National University in June 2003. In early 2005 she spent three months as a visiting student at the IBM Tokyo Research Laboratory. Her current research interests are in optimization algorithm development with applications to control design and related problems.
Robert Orsi received his Ph.D. in 2000 from the University of Melbourne, Australia. Since then he has held positions with Voyan Technology, the Australian National University, and National ICT Australia. Currently, since July 2004, he is the recipient of an Australian Research Council Australian Postdoctoral Fellowship, which he has taken up at the Australian National University. His current research interests include optimization methods for engineering problems, convex optimization, rank constrained LMIs, and inverse eigenvalue problems.
Brief paper
Generalized pole placement via static output feedback: A methodology based on projections夡 Kaiyang Yang a,b,∗ , Robert Orsi a a Research School of Information Sciences and Engineering, The Australian National University, Canberra ACT 0200, Australia b National ICT Australia Limited, Locked Bag 8001, Canberra ACT 2601, Australia
Received 30 September 2005; received in revised form 14 March 2006; accepted 29 June 2006 Available online 30 August 2006
Abstract This paper presents an algorithm for solving static output feedback pole placement problems of the following rather general form: given n subsets of the complex plane, find a static output feedback that places in each of these subsets a pole of the closed-loop system. The algorithm presented is iterative in nature and is based on alternating projection ideas. Each iteration of the algorithm involves a Schur matrix decomposition, a standard least-squares problem and a combinatorial least-squares problem. While the algorithm is not guaranteed to always find a solution, computational results are presented demonstrating the effectiveness of the algorithm. 䉷 2006 Elsevier Ltd. All rights reserved. Keywords: Static output feedback; Pole placement; Stabilization; Alternating projections
1. Introduction
In this paper we will actually consider the following generalized static output feedback pole placement problem.
There has been a great deal of research done on the problems of pole placement and stabilization via static output feedback. An overview of theoretical results, existing algorithms and historical developments can be found in Byrnes (1989), Rosenthal and Willems (1998), Syrmos, Abdallah, Dorato, and Grigoriadis (1997), de Oliveira and Geromel (1997), Rosenthal and Sottile (1998) and Eremenko and Gabrielov (2002). Our main interest here is in algorithms, and in this regard, for pole placement, the survey paper Rosenthal and Willems (1998) states that existing sufficiency conditions are mainly theoretical in nature and that there are no good numerical algorithms available in many cases when a problem is known to be solvable. Despite the great deal of work that has been done in this area, new algorithms for these important problems are still of interest.
Problem 1. Given A ∈ Rn×n , B ∈ Rn×m , C ∈ Rp×n and closed subsets C1 , . . . , Cn ⊂ C, find K ∈ Rm×p such that i (A + BKC) ∈ Ci for i = 1, . . . , n.
夡 A preliminary version of this paper was presented at the 16th IFAC World Congress. This paper was recommended for publication in revised form by Associate Editor Yasumasa Fujisaki under the direction of Editor Roberto Tempo. ∗ Corresponding author. Research School of Information Sciences and Engineering, The Australian National University, Canberra ACT 0200, Australia. E-mail addresses: [email protected] (K. Yang), [email protected] (R. Orsi).
0005-1098/$ - see front matter 䉷 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.automatica.2006.06.021
Here i (A + BKC) denotes the ith eigenvalue of A + BKC. Problem 1 encompasses many types of pole placement problems. Here are some examples: (1) Classical pole placement: Ci = {ci },
ci ∈ C.
(2) Stabilization type problems for continuous time and discrete time systems: C1 = · · · = Cn = {z ∈ C | Re(z) − },
>0
and C1 = · · · = Cn = {z ∈ C | |z| }, respectively. (3) Relaxed classical pole placement: Ci = {z ∈ C | |z − ci | ri }.
0 < < 1,
(1)
2144
K. Yang, R. Orsi / Automatica 42 (2006) 2143 – 2150
Fig. 1. A hybrid problem.
Here each region Ci is a disk centered at ci ∈ C with radius ri 0. (4) Hybrid problems: For example, problems of the type shown in Fig. 1. Here c ∈ C and the aim is to place a pair of poles at c and c, ¯ and to place the remaining poles in the truncated cone C: C1 = {c},
C2 = {c}, ¯
C3 = · · · = Cn = C.
As far as the authors are aware, pole placement in the generality presented in Problem 1 has not previously been considered. The most closely related results from the existing literature are as follows. Early work in Gutman and Jury (1981) considers pole placement in a single region specified by polynomials. While a Lyapunov-type necessary and sufficient condition is given for a matrix to have its eigenvalues in such a region, this condition is polynomial in the matrix in question and hence not readily amenable to controller design. In Chilali and Gahinet (1996), LMI conditions are presented that are sufficient though not necessary for pole placement in various convex regions. These results cover state feedback and fullorder dynamic output feedback but not static output feedback. Another LMI-based approach, again for a single, though this time, possibly disconnected region, is considered in Bosche, Bachelier, and Mehdi (2004). A method for placing poles in distinct convex regions (each region is specified using linear programming constraints or second-order cone constraints) is given in Hassibi, How, and Boyd (1999); however, the method is based on eigenvalue perturbation results and hence appears largely limited to cases where the open-loop poles are already quite close to the desired poles. In Satoh, Okubo, and Sugimoto (2003), pole placement in distinct convex regions (each region is a disk or a half plane) is achieved via a rank constrained LMI approach though the results are only for state feedback. This paper presents an algorithm for Problem 1. The approach employed here is quite different to each of the approaches mentioned above. Problem 1 is shown to be equivalent to finding a point in the intersection of two particular sets, one of which is a simple convex set and the other a rather complicated nonconvex set. The algorithm is iterative in nature and is based on an alternating projection like scheme between these two sets. Each iteration of the algorithm involves a Schur matrix decomposition and a standard least-squares problem. If the Ci ’s are not all equal, each iteration also requires a combinatorial least-squares matching step.
Alternating projection type ideas have been employed previously for output feedback stabilization, see in particular Grigoriadis and Skelton (1996), Grigoriadis and Beran (1999) and Orsi, Helmke, and Moore (2006).1 A distinguishing feature of our algorithm is that, unlike these methods, our algorithm does not involve LMIs. (A further technical difference is that rather than solving feasibility problems that involve symmetric matrices, the problem solved by the algorithm is a feasibility problem involving nonsymmetric matrices.) The algorithm can be applied to problems with rather general choices for the Ci regions. In fact the only formal requirement is the following, which we state in the form of an assumption. Assumption 1. It is possible to calculate projections onto each of the Ci ’s: given z ∈ C it is possible to find zi ∈ Ci such that |z − zi | |z − c| for all c ∈ Ci . In particular, the Ci ’s must be closed sets though they need not be convex or even connected. For given A, B, C and Ci ’s, Problem 1 may or may not have a solution. Indeed, one would expect that determining whether a particular instance of Problem 1 is solvable is in general difficult. For example, the problem of determining whether the classical pole placement problem is solvable for particular A, B, C and desired poles has recently been shown to be NPhard (Fu, 2004). Given the difficulty of Problem 1, an efficient (i.e., polynomial time) algorithm that is able to correctly solve all instances of the problem cannot be expected, and while the algorithm presented here is often quite effective in practice, it is not guaranteed to find a solution even if a solution exists. The structure of the paper is as follows. The remainder of this section lists some notation which is used in the rest of the paper. Projections play a key part in the algorithm and Section 2 contains general properties of projections and recalls how alternating projections can be used to find a point in the intersection of a finite number of closed (convex) sets. Section 3 presents the solution methodology. To motivate the solution methodology we first restrict our attention to systems that have a symmetric state space representation (A = AT , C = B T ) and present an algorithm for this easier class of problems. The general problem is then considered. Section 4 contains computational results of applying the algorithm to various instances of Problem 1. The paper ends with some concluding remarks. Notation: Sr is the set of real symmetric r × r matrices. diag(v) for v ∈ Cr denotes the r × r diagonal matrix whose ith diagonal term is vi . For Z ∈ Cr×s , vec(Z) ∈ Crs consists of the columns of Z stacked below each other. Y ⊗ Z denotes the Kronecker product of Y and Z. For Z ∈ Cr×s , Re(Z) ∈ Rr×s and Im(Z) ∈ Rr×s denote, respectively, the real and imaginary parts of Z. 1 They were first used in control design to solve certain convex problems, see Grigoriadis and Skelton (1994).
K. Yang, R. Orsi / Automatica 42 (2006) 2143 – 2150
2. Projections Projections play a key part in the algorithm. This section contains some general properties of projections that are used throughout the paper. Let x be an element in a Hilbert space H and let D be a closed (possibly nonconvex) subset of H. Any d0 ∈ D such that x − d0 x − d for all d ∈ D will be called a projection of x onto D. In the cases of interest here, namely that H is a finite dimensional Hilbert space, there is always at least one such point for each x. If D is convex as well as closed then each x has exactly one such minimum distance point. A function PD : H → H will be called a projection operator (for D) if for each x ∈ H , PD (x) is a projection of x onto D. In this paper the problems we are interested in solving are feasibility problems of the following general form. Problem 2. Given closed sets D1 , . . . , DN in a finitedimensional Hilbert space H, find a point in the intersection N i=1 Di (assuming the intersection is non-empty). (In fact we will only be interested in the case N = 2.) When all the Di ’s in Problem 2 are convex, a point in the intersection of the Di ’s can be found via alternating projections. Theorem 3. Let D1 , . . . , DN be closed convex sets in a finite dimensional Hilbert space H. Given constants t1 , . . . , tN in the interval (0, 2), for i = 1, . . . , N, define the operators Ri = (1 − ti )I d + ti PDi . If N i=1 Di is non-empty, then starting from an arbitrary initial value, the sequence defined by xi+1 = Ri (mod N)+1 (xi ) converges to an element in N D . i i=1
projections onto these sets are computable, a solution method is to alternately project onto D1 and D2 . Theorem 4 ensures that the distance xn − yn is nonincreasing with n. While the convergence of xn −yn to 0 is not guaranteed, if this quantity does converge to 0 then a solution to the problem will have been found. 3. Methodology 3.1. The symmetric problem Systems with a symmetric state space realization, that is, systems with state space matrices satisfying A = AT , C = B T , occur in various contexts, for example RC-networks. In order to motivate our solution methodology, we first consider the following special case of Problem 1 for symmetric systems.2 Problem 5. Given A ∈ Sn , B ∈ Rn×m , and closed subsets C1 , . . . , Cn ⊂ R, find K ∈ Sm such that i (A + BKB T ) ∈ Ci for i = 1, . . . , n. Note that in Problem 5 the static output feedback matrix K is required to be symmetric and that the Ci ’s are assumed to be subsets of the real numbers. This latter assumption is natural as given any symmetric K, A + BKB T is also symmetric and hence has real eigenvalues. We assume A, B and C1 , . . . , Cn are given and fixed. Let L = {Z ∈ Sn | Z = A + BKB T for some K ∈ Sm } and let M denote the set of symmetric matrices with eigenvalues in the specified regions C1 , . . . , Cn ,
Proof. See Gubin, Polyak, and Raik (1967). An alternate proof can also be found in Youla and Webb (1982).
M = {Z ∈ Sn | i (Z) ∈ Ci , i = 1, . . . , n}.
Here Id denotes the identity operator on H; I d(x) = x for all x ∈ H . Each Ri of the theorem is a projection operator onto Di if and only if ti = 1. If ti = 1 then Ri (x) is referred to as a relaxed projection of x. Theorem 3 is a generalization of the classical method of alternating projections which corresponds to the case t1 = · · · = tN = 1. We will see later how the freedom to choose the ti ’s not equal to 1 can be useful. When one or more Di ’s are nonconvex, Theorem 3 no longer applies and starting the algorithm of Theorem 3 from certain initial values may result in a sequence of points that does not converge to a solution of the problem (Combettes & Trussell, 1990). If, however, there are only two sets, then the following distance reduction property always holds (Orsi & Yang, 2006).
The symmetric problem can be stated as follows:
Theorem 4. Let D1 and D2 be closed, nonempty, possibly nonconvex sets in a finite dimensional Hilbert space H, and let PD1 and PD2 be associated projection operators. For any initial value y0 ∈ H and for all n 0, let xn+1 = PD1 (yn ) and yn+1 =PD2 (xn+1 ). Then for all n 1, xn+1 −yn+1 xn+1 − yn xn − yn . Suppose one is interested in solving Problem 2 in the case of two sets, D1 and D2 , when one or both sets are nonconvex. If
2145
Find X ∈ L ∩ M. We now show that projections onto both sets L and M can be calculated and hence that alternating projections can be employed as a solution method. From now on, Sn will be regarded as a Hilbert space with inner product Y, Z = tr(Y Z) and associated norm Z = Z, Z 1/2 (the Frobenius norm). For i = 1, . . . , n, PCi will denote a projection operator for Ci . The set L is an affine subspace and hence is convex. Projection of X ∈ Sn onto L involves solving a standard leastsquares problem. The following result is fairly straightforward to establish. Lemma 6. The projection of X ∈ Sn onto L is given by PL (X) = A + BKB T where K ∈ Sm is a minimizer of (B ⊗ B)vec(K) − vec(X − A)2 . ( · 2 denotes the standard vector 2 norm.) 2 See Mahony and Helmke (1995) for more regarding the classical pole placement problem for symmetric systems and a rather surprising result regarding arbitrary pole placement for such systems.
2146
K. Yang, R. Orsi / Automatica 42 (2006) 2143 – 2150
While the set M is not convex, it is still possible to calculate projections onto M. How to do this is shown in Theorem 8. The result is based on the following result of Hoffman and Wielandt. Lemma 7. Suppose Y, Z ∈ Sn have eigenvalue-eigenvector decompositions Y = V DV T , D = diag(Y1 , . . . , Yn ), and Z = Z n×n W EW T , E = diag(Z are or1 , . . . , n ), where V , W ∈ R Y Y Z Z thogonal and 1 · · · n and 1 · · · n . Then D − EY − Z,
(2)
where · denotes the Frobenius norm.
Theorem 8. Given Y ∈ Sn , let Y = V DV T be an eigenvalueeigenvector decomposition of Y with D = diag(1 , . . . , n ). Let be a permutation of {1, . . . , n}such that, amongst all possible permutations, it minimizes nk=1 |k − PC(k) (k )|2 . ˆ T where Dˆ = diag(PC (1 ), . . . , Define PM (V , D) = V DV (1) PC(n) (n )). Then PM (V , D) is a best approximant in M to Y in the Frobenius norm. Proof. Let Y be as in the theorem statement. As PM (V , D) ∈ M, it remains to show for all Z ∈ M.
(3)
Without loss of generality, suppose the eigenvalues of Y are ordered, i.e., 1 · · · n . Similarly, for Z ∈ M, let Z = W EW T be an eigenvalue-eigenvector decomposition with E = Z Z Z diag(Z 1 , . . . , n ) and 1 · · · n . As the Frobenius norm is orthogonally invariant, ˆ Y − PM (V , D) = D − D.
(4)
Let be a permutation of {1, . . . , n} such that Z k ∈ C(k) , k = 1, . . . , n. (Such a permutation exists as Z ∈ M and hence must have an eigenvalue in each of the Ci ’s.) Define Eˆ = diag(PC(1) (1 ), . . . , PC(n) (n )). It follows from the definition of Dˆ that ˆ D − E. ˆ D − D
(5)
As for each k, |k − PC(k) (k )||k − Z k |, it also follows that ˆ D − ED − E.
L = {Z ∈ Rn×n | Z = A + BKC for some K ∈ Rm×p }, and re-define M to be the set of complex matrices with eigenvalues in the specified regions C1 , . . . , Cn , M = {Z ∈ Cn×n | i (Z) ∈ Ci , i = 1, . . . , n}. Problem 1 can now be stated as
Proof. See for example Horn and Johnson (1985, Corollary 6.3.8).
Y − PM (V , D) Y − Z
From now on Cn×n will be regarded as a Hilbert space with inner product Y, Z = tr(Y Z ∗ ). In this subsection we re-define L to be the set of all possible closed-loop system matrices,
(6)
Combining (4)–(6) and inequality (2) from Lemma 7 gives the inequality in (3) and the proof is complete. Note that to calculate PM (V , D) we keep the original orthogonal matrix V and simply modify the diagonal matrix D to ˆ The fact that V remains unchanged motivates our solution D. method for the general nonsymmetric case. 3.2. The general nonsymmetric problem Consider again Problem 1. Throughout this subsection it is assumed A, B, C and C1 , . . . , Cn are given and fixed.
Find X ∈ L ∩ M. A solution strategy to solve Problem 1 would be to employ an alternating projection scheme, alternatively projecting between L and M. A difficulty occurs in trying to do this. While L is convex and M is nonconvex, just as for the symmetric problem, M in this case is a substantially more complicated set and how to calculate projections onto this set is a difficult unsolved problem. As will be verified by the experiments, an alternating projection like scheme can still be quite successful if for M a suitable substitute mapping is used in place of an actual projection operator. Before proceeding, recall Schur’s result (Horn & Johnson, 1985). Theorem 9. Given Z ∈ Cn×n with eigenvalues 1 , . . . , n in any prescribed order, there is a unitary matrix V ∈ Cn×n and an upper triangular matrix T ∈ Cn×n such that Z = V T V ∗ , and Tkk = k for k = 1, . . . , n. The following map is proposed as a substitute for a projection map onto M. Though it is not a true projection map, the notation PM will still be used. Our choice for PM is motivated by the projection operator for M for the symmetric problem and the idea for using such a substitute mapping originates in Orsi and Yang (2006) (see also Orsi (2006)). We stress that as PM is not a true projection map, unlike for the symmetric problem, the distance reduction property of Theorem 4 may not hold. While PM has various desirable properties (see below), the theoretical convergence properties of the algorithm are currently unclear, providing interesting questions for future research. In what follows, Assumption 1 is assumed to hold. Definition 10. Suppose V ∈ Cn×n is unitary and T ∈ Cn×n is upper triangular. Let be a permutation of {1, . . . , n} such that amongst all possible permutations, it minimizes n
|Tkk − PC(k) (Tkk )|2 .
(7)
k=1
Define PM (V , T ) = V Tˆ V ∗ where Tˆ is upper triangular and given by PC(k) (Tkk ) if k = l, Tˆkl = Tkl otherwise.
K. Yang, R. Orsi / Automatica 42 (2006) 2143 – 2150
Note that PM maps into the set M. Note also that, just as for the symmetric problem, finding involves solving a combinatorial least-squares problem. (This will be discussed further later in the section.) In order to apply PM to Z ∈ Cn×n , a Schur decomposition of Z must first be found. A given Z may have a nonunique Schur decomposition and Z = V1 T1 V1∗ = V2 T2 V2∗ does not necessarily imply PM (V1 , T1 ) = PM (V2 , T2 ). Hence, PM may give different points for different Schur decompositions of the same matrix. (Though it was not discussed at the time, similar comments apply to the PM mapping for the symmetric case.) This is not so important as different Schur decompositions lead to points in M of equal distance to the original matrix. Also note that of all the points in M that have a Schur decomposition with the same V matrix, PM (V , T ) is closest (or at least equal closest) to the original point Z = V T V ∗ . For a proof of these last two claims, see Orsi and Yang (2006) or Orsi (2006). As for the symmetric problem, projection of X ∈ Cn×n onto L involves solving a standard least-squares problem. Lemma 11. The projection of X ∈ Cn×n onto L is given by PL (X) = A + BKC where K ∈ Rm×p is a minimizer of (C T ⊗ B)vec(K) − vec(Re(X) − A)2 . Here is our algorithm for Problem 1. Algorithm. Problem Data: A ∈ Rn×n , B ∈ Rn×m , C ∈ Rp×n and C1 , . . . , Cn ⊂ C. Initialization: Choose a randomly generated matrix Y ∈ Rn×n . For example, draw each entry of Y from a normal distribution of zero mean and variance 1. repeat 1. X := PL (Y ). 2. Calculate a Schur decomposition of X: X = V T V ∗ . 3. Y := PM (V , T ). until X − Y < . see Definition 10, Note that as Y = PM (V , T ) = V Tˆ V ∗ , X −Y =V T V ∗ −V Tˆ V ∗ =T − Tˆ =( |Tkk − Tˆkk |2 )1/2 . As the Tkk ’s are the eigenvalues of X ∈ L and the Tˆkk ’s are the eigenvalues of Y ∈ M, the algorithm stops when X (which equals A+BKC for some K) has eigenvalues sufficiently close to those of a matrix that satisfies the pole placement constraints. In particular, each eigenvalue of such an X cannot violate the pole placement constraints by more than . As mentioned previously, PM involves finding a permutation that minimizes (7). The first step in solving this combinatorial least-squares problem is calculating the squared distance of each Tkk to each subset Cl and placing these values in an n × n cost matrix D: Dkl := |Tkk − PCl (Tkk )|2 .
(8)
The problem is now equivalent to finding a permutation such that Dk (k) is minimal.
2147
The problem of finding a minimizing given a cost matrix D is a linear assignment problem which can be solved in O(n3 ) time using the so-called Hungarian method, see Martello and Toth (1987) for details. Note that if the Ci ’s are not all distinct, for example this occurs for stabilization problems and what here have been called hybrid problems, the complexity of the matching problem is reduced. In fact for stabilization problems, all the Ci ’s are the same and no matching step is required. For the hybrid problem shown in Fig. 1, it is only necessary to check n(n − 1) possibilities corresponding to which two Tkk ’s are matched to c and c. ¯ Hence for this hybrid problem, the direct approach is faster than using the Hungarian method. Also note that, given a cost matrix D, an alternative to the Hungarian method is the following faster suboptimal matching strategy. Find the (or a) smallest entry in D, say Dk˜ l˜. Match Tk˜ k˜ with Cl˜ and cross out row k˜ and column l˜ of D. Now only consider the uncrossed out entries in D and repeat, until all n matches have been made. This method does not always find the optimal matching though it can often be quite an effective substitute for the Hungarian method. It will be termed suboptimal matching. Surprisingly, as will be shown in the next section, by using suboptimal matching in the algorithm it was possible to solve a particular problem which was not solvable using the Hungarian method. 4. Computational results This section contains computational results of applying the algorithm to various instances of Problem 1. We include results for classical pole placement, continuous time stabilization, and a hybrid problem. (Additional results can be found in Yang & Orsi (2005) and Yang, Orsi, & Moore (2004).) The algorithms for each problem were implemented in Matlab 6.5 and all results were obtained using a 3.06 GHz Pentium 4 machine. Throughout this section a randomly generated matrix will be a matrix whose entries are drawn from a normal distribution of zero mean and variance 1. 4.1. Classical pole placement: random problems This subsection contains results for some randomly generated classical pole placement problems. A 1000 problems with n=6, m = 4 and p = 3 were created. Each problem was created as follows. A, B, C and K˜ were generated randomly. A scalar multiple of the identity was added to A to ensure the stability ˜ was equal to = 0.1. The desired poles degree of A + B KC ˜ were taken to be the poles of A + B KC. An attempt was made to solve each problem using up to 10 different initial conditions and a maximum of a 1000 iterations per initial condition. With the termination parameter set to = 10−3 , the overall success rate was 91%. (The success rate increases to 95% if up to 20 initial conditions are used.) For the problems that were solved, the average number of iterations taken was 1.8 × 103 and the average time taken was 1.5 CPU seconds.
2148
K. Yang, R. Orsi / Automatica 42 (2006) 2143 – 2150
Table 1 A comparison of performance for different n, m and p. i denotes the number of iterations and listed are the number of problems that converged within different iteration ranges. ‘NC’ denotes the number of problems that did not converge within 1000 iterations. T denotes the average convergence time in CPU seconds n m p
3 1 1
6 1 1
9 1 1
3 2 1
6 4 3
9 5 4
1 i 10 10 < i 100 100 < i 1000 NC
996 3 1 0
990 7 3 0
991 8 1 0
993 3 4 0
989 11 0 0
982 17 1 0
T
0.0007
0.0007
0.001
0.0006
0.001
0.002
4.2. Classical pole placement: a particular problem from the literature The following problem is taken from Sridhar and Lindorff (1973) and is of interest as the set of desired poles overlaps with the set of open-loop poles. The system matrices are the following: ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ 1 0 0 0 1 0 1 0 T 0 ⎥ ⎢0 2 0 ⎢0 1⎥ ⎢1 0⎥ A=⎣ ⎦, B = ⎣ ⎦, C = ⎣ ⎦ . 0 0 −3 0 1 0 0 1 1 1 0 1 0 0 0 −4 The aim is to place the closed-loop poles at {−1, −2, −3, −5}. While initial attempts to solve this problem failed, the problem was solved by using suboptimal matching and replacing Step 3 of the algorithm with the relaxed projection ‘Y := (1 − t)X + tP M (V , T )’, t ∈ (0, 2) constant. Strictly speaking, this is not a relaxed projection in the true sense as PM (V , T ) may not be a projection of X = V T V ∗ , however, the idea is clearly the same. Solutions were successfully found by taking t close to 0; t =0.1, 0.2 and 0.3 can all be used to successfully find a solution. Likelihood of success increases and speed of convergence decreases the closer t is to 0. With t = 0.3 and = 10−3 , solutions can typically be found in about 1.2 × 104 iterations and about 4.7 CPU seconds. Note: When employing relaxed projections, the loop termination criterion of the algorithm should be replaced by ‘until X − PM (V , T ) < ’ as now X − Y = X − ((1 − t)X + tP M (V , T )) = tX − PM (V , T ). 4.3. Continuous time stabilization: random problems In the comparison paper de Oliveira and Geromel (1997), a number of methods for continuous time stabilization via static output feedback were compared. In this subsection we repeat the same numerical experiments of de Oliveira and Geromel (1997) using our algorithm, to see how it compares. One thousand problems were generated for each of a number of different choices of system dimensions n, m and p. In each case, A, B, C and K˜ were generated randomly. A scalar multiple of the identity was added to A to ensure it had a stability degree ˜ All stable A’s were of one. A was then replaced by A − B KC. discarded.
In applying our algorithm, the Ci ’s were chosen as in (1) with = 0.1. As closed-loop systems were only required to be stable, the algorithm was terminated as soon as all poles had real part less than zero. An attempt was made to solve each problem using up to 10 different initial conditions and a maximum of a 100 iterations per initial condition. Results are given in Table 1. As can be seen, all problems were solved. Most problems were solved in 10 iterations or less and average solution times were very small. These results are as good as the results for the best two algorithms tested in de Oliveira and Geromel (1997). Other values for , such as 0.2, 0.5, 1 and 2, were also tried and produced just as good results, indicating a certain robustness of the algorithm with respect to this parameter. While the above results are very encouraging, the algorithm may not work as well for all problems. For example, in Yang et al. (2004) the method of problem generation is different and the algorithm does not perform as well. 4.4. A hybrid problem So far we have presented results for two traditional classes of problems. In this subsection we demonstrate the generality of the algorithm by considering a less standard problem, namely a hybrid problem of the type shown in Fig. 1. The problem parameters are c = −0.5 + i3, C = {z ∈ C | Re(z) − 2 and |Im(z)| |Re(z)|}, n = 13, m = 3 and p = 5. To ensure solvability, B, C and K˜ were ran˜ with domly generated and A set to A = V˜ T˜ V˜ T − B KC ˜ ˜ V a random orthogonal matrix and T a real block upper triangular matrix with a spectrum satisfying the constraints. (T˜ was assigned the spectrum {−0.5 ± i3, −2, −2 ± i, −2.3, −2.5, −3 ± i3, −3.5 ± i3.1, −4 ± i4} and was created by choosing appropriate 1 × 1 and 2 × 2 blocks for its diagonal. The remaining upper triangular entries of T˜ were chosen randomly.) With = 10−3 , 64% of initial conditions tested resulted in a solution of this problem within 5000 iterations. For the initial conditions that lead to convergence, the average number of iterations taken was 8.6 × 102 and the average time taken was 3.1 CPU seconds. The closed-loop poles corresponding to a particular solution are shown in Fig. 2.
K. Yang, R. Orsi / Automatica 42 (2006) 2143 – 2150 5 4 3 2 1 0 -1 -2 -3 -4 -5
-5
-4.5
-4
-3.5
-3
-2.5
-2
-1.5
-1
-0.5
0
Fig. 2. The closed-loop poles corresponding to a solution for the considered hybrid problem.
Note: The least-squares matching steps were done directly rather than using the Hungarian method. 5. Conclusion In this paper a new methodology for solving a broad class of output feedback pole placement problems has been presented. While the methodology is not guaranteed to find a solution, numerical experiments presented in the paper demonstrate it can be quite effective in practice. A particular strength of the algorithm is that, in addition to being able to solve classical pole placement and stabilization problems, it can also be used to solve less standard pole placement problems. Acknowledgments The authors thank Iven Mareels and Robert Mahony for their helpful guidance to various parts of the pole placement literature. The second author acknowledges the support of the Australian Research Council through Grant DP0450539. National ICT Australia is funded through the Australian Government’s Backing Australia’s Ability initiative, in part through the Australian Research Council. References Bosche, J., Bachelier, O., & Mehdi, D. (2004). Robust pole placement by static output feedback. In Proceedings of the 43rd IEEE conference on decision and control (pp. 869–874). Paradise Island, Bahamas. Byrnes, C. I. (1989). Pole assignment by output feedback. In Nijmeijer, & Schumacher (Eds.), Three decades of mathematical system theory, Lecture notes in control and information science Vol. 135, (pp. 31–78). Berlin: Springer. Chilali, M., & Gahinet, P. (1996). H∞ design with pole placement constraints: An LMI approach. IEEE Transactions on Automatic Control, 41(3), 358–367. Combettes, P. L., & Trussell, H. J. (1990). Method of successive projections for finding a common point of sets in metric spaces. Journal of Optimization Theory and Applications, 67, 487–507.
2149
de Oliveira, M. C., & Geromel, J. C. (1997). Numerical comparison of output feedback design methods. In Proceedings of American control conference (pp. 72–76). New Mexico, USA. Eremenko, A., & Gabrielov, A. (2002). Pole placement by static output feedback for generic linear systems. SIAM Journal of Control Optimization, 41(1), 303–312. Fu, M. (2004). Pole placement via static output feedback is NP-hard. IEEE Transactions on Automatic Control, 49(5), 855–857. Grigoriadis, K. M., & Beran, E. B. (1999). Alternating projection algorithms for linear matrix inequalities problems with rank constraints. In L. El Ghaoui, & S.-I. Niculescu (Eds.), Advances on linear matrix inequality methods in control (pp. 251–267). Philadelphia, PA: SIAM. Grigoriadis, K. M., & Skelton, R. E. (1994). Alternating convex projection methods for covariance control design. International Journal of Control, 60(6), 1083–1106. Grigoriadis, K. M., & Skelton, R. E. (1996). Low-order control design for LMI problems using alternating projections. Automatica, 32(8), 1117–1125. Gubin, L. G., Polyak, B. T., & Raik, E. V. (1967). The method of projections for finding the common point of convex sets. USSR Computational Mathematics and Mathematical Physics, 7, 1–24. Gutman, S., & Jury, E. (1981). A general theory for matrix root-clustering in subregions of the complex plane. IEEE Transactions on Automatic Control, 26(4), 853–863. Hassibi, A., How, J., & Boyd, S. (1999). A path-following method for solving BMI problems in control. In Proceedings of the American control conference (Vol. 2) (pp. 1385–1389), San Diego, USA. Horn, R. A., & Johnson, C. R. (1985). Matrix analysis. Cambridge: Cambridge University Press. Mahony, R. E., & Helmke, U. (1995). System assignment and pole placement for symmetric realisations. Journal of Mathematical Systems, Estimation and Control, 5(2), 1–32. Martello, S., & Toth, P. (1987). Linear assignment problems. Annals of Discrete Mathematics, 31, 259–282. Orsi, R. (2006). Numerical methods for solving inverse eigenvalue problems for nonnegative matrices. SIAM Journal on Matrix Analysis and Applications, 28(1), 190–212. Orsi, R., Helmke, U., & Moore, J. B. (2006). A Newton-like method for solving rank constrained linear matrix inequalities. Automatica, 42 (11). Orsi, R., & Yang, K. (2006). Numerical methods for solving inverse eigenvalue problems for nonnegative matrices. In Proceedings of the 17th international symposium on mathematical theory of networks and systems (pp. 2284–2290). Kyoto, Japan. Rosenthal, J., & Sottile, F. (1998). Some remarks on real and complex output feedback. Systems and Control Letters, 33, 73–80. Rosenthal, J., & Willems, J. C. (1998). Open problems in the area of pole placement. In V.D. Blondel, E.D. Sontag, M. Vidyasager, & J.C. Willems (Eds.), Open problems in mathematical systems and control theory (pp. 181–191). Berlin: Springer. Satoh, A., Okubo, J., & Sugimoto, K. (2003). Transient response shaping in H∞ control by eigenstructure assignment to convex regions. In Proceedings of the 42nd IEEE conference on decision and control (pp. 2288–2293). Maui, USA. Sridhar, B., & Lindorff, D. P. (1973). Pole placement with constraint gain output feedback. International Journal of Control, 18(5), 993–1003. Syrmos, V. L., Abdallah, C. T., Dorato, P., & Grigoriadis, K. (1997). Static output feedack—A survey. Automatica, 33(2), 125–137. Yang, K., & Orsi, R. (2005). Pole placement via output feedback: A methodology based on projections. In Proceedings of the 16th IFAC world congress, Prague, Czech Republic. Yang, K., Orsi, R., & Moore, J. B. (2004). A projective algorithm for static output feedback stabilization. In Proceedings of the second IFAC symposium on system, structure and control (pp. 263–268). Oaxaca, Mexico. Youla, D. C., & Webb, H. (1982). Image restoration by the method of convex projections: Part 1—Theory. IEEE Transactions on Medical Imaging, MI1(2), 81–94.
2150
K. Yang, R. Orsi / Automatica 42 (2006) 2143 – 2150 Kaiyang Yang was born in Shenyang, PR China, in 1979. She received her B.E. degree in Control Science and Engineering from Zhejiang University, PR China, in 2002. From July 2002 to June 2003, she worked for Jinboyuan System Integration as a system engineer. Ms. Yang started her Ph.D. studies at the Australian National University in June 2003. In early 2005 she spent three months as a visiting student at the IBM Tokyo Research Laboratory. Her current research interests are in optimization algorithm development with applications to control design and related problems.
Robert Orsi received his Ph.D. in 2000 from the University of Melbourne, Australia. Since then he has held positions with Voyan Technology, the Australian National University, and National ICT Australia. Currently, since July 2004, he is the recipient of an Australian Research Council Australian Postdoctoral Fellowship, which he has taken up at the Australian National University. His current research interests include optimization methods for engineering problems, convex optimization, rank constrained LMIs, and inverse eigenvalue problems.