# On computing stabilizability radii of linear time-invariant continuous systems.

1. Introduction. Consider the linear time-invariant continuous dynamical system

(1.1) x(t) = Ax(t) + Bu(t), t [greater than or equal to] 0,

where A [member of] [C.sup.nxn], B [member of] [C.sup.nxm], with state vector x(t) [member of] [C.sup.n] and control vector u(t) [member of] [C.sup.m] for all t [greater than or equal to] 0. The system (1.1) is called stabilizable if there exists a feedback u(t) = Fx(t), where F is a fixed matrix, making the closed-loop system x(t) = (A + BF)x(t) asymptotically stable. It is well-known that the system (1.1) is stabilizable if and only if the condition rank [A-[lambda]I B]= n holds for all [lambda] [member of] [C.sub.+] := {[lambda] [member of] C : R[lambda] [greater than or equal to] 0}; see, for example, .

One of the most effective and flexible approaches towards problems of robustness of stabilizability is based on the concept of the stabilizability radius--the norm of the smallest perturbation that makes a given system unstabilizable, introduced in [5, 6] as

[[tau].sub.0] (A,B) := inf {[parallel][[[DELTA].sub.A], [[DELTA].sub.B]][parallel] : the perturbed system

x = (A + [[DELTA].sub.A])x + (B + [[DELTA].sub.B])u is unstabilizable} ,

where [parallel]*[parallel] denotes the spectral norm. This definition is inspired from the definition of the controllability radius in . The stabilizability radius is expressed in [5, 6] as

(1.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [[sigma].sub.min](*) denotes the smallest singular value of its matrix argument. Recently in , we considered stabilizability radii when only one of the system matrices, A or B, is perturbed

[[tau].sub.1](A) := inf {[parallel] [DELTA].sub.A] [parallel] : the perturbed system x = (A + [[DELTA].sub.A])x + Bu is unstabilizable} ,

[[tau].sub.2] (B) := inf {[parallel] [[DELTA].sub.B] [parallel] : the perturbed system x = Ax + (B + [[DELTA].sub.B])u is unstabilizable},

and the formulas for these values are given by

(1.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(1.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where null (*) is the matrix whose columns form an orthonormal basis of the null space of its matrix argument.

Problem (1.4) can be computed easily. But problems (1.2) and (1.3) are non-smooth optimization problems in two real variables [alpha] and [beta], the real and imaginary parts of [lambda]. Moreover, the objective functions [[sigma].sub.min] ([A - [lambda]I B]) or [[sigma].sub.min] ([null (B*)]*(A - [lanbda]I)) are not convex and may have many local minima, so standard optimization methods, which usually are only guaranteed to converge to a local minimum, will not yield reliable results in general.

To the best of our knowledge, the problem of computing stabilizability radii has not been studied in the literature even though several algorithms have been designed to compute controllability radii; see [3, 4, 7, 8, 9, 11, 13]. In the papers [4, 7, 11], two-dimensional grid techniques are used that are too costly for high accuracy. Gu  proposed a bisection method which can correctly estimate [[tau].sub.0] (A, B) within a factor of two in polynomial time in n. Burke et al.  suggested a trisection variant to retrieve the distance to uncontrollability to any desired accuracy with O([n.sup.6]) complexity. Gu et al.  and Mengi  reduced the average running time to O([n.sup.4]) by employing inverse iterations and the shift-and-invert preconditioned Arnoldi method. It is the main purpose of this paper to describe a numerical method for computing both problems (1.2) and (1.3).

Indeed, we consider the following non-convex and non-smooth general optimization problem

(1.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where P, Q are some given matrices in [C.sup.pxq] such that p [less than or equal to] q and rank(Q) = p.

We can check that if P := [A B], Q := [I 0], then problem (1.5) reduces to problem (1.2), and if P := [null(B*)]*A, Q := [null (B*)]*, then problem (1.5) reduces to problem (1.3). In this paper, based on the idea of the trisection algorithm introduced in , we present a method to solve (1.5).

The structure of this paper is as follows. In the next section, we give a modified version of Gu's result [8, Theorem 3.1] that is applicable to (1.5). Then, we apply the obtained results to state a method for solving this problem in Section 3. Finally, the reliability of the algorithm is demonstrated by numerical examples in Section 4.

2. Modified version of Gu's theorem. The methods for computing the controllability radius in [8, 9, 13] are based on a simultaneous comparison of two bounds [[delta].sub.1] > [[delta].sub.2] with [[tau].sub.0] (A, B), i.e., one of the following inequalities is verified

[[tau].sub.0](A,B) [less than or equal to] [[delta].sub.1] or [[tau].sub.0](A,B) > [[delta].sub.2].

This so-called Gu's test based on [8, Theorem 3.1] returns some information about only one of the inequalities even if both of them may be satisfied.

It is remarkable that Gu's test in  cannot be applied to solve the general problem (1.5). The search space of (1.5) is just the closed right half of the complex plane, which naturally leads to the idea of a vertical search as in the following theorem.

THEOREM 2.1. Assume that [delta] > [tau](P, Q). Given a number [eta] [member of] (0, 2([delta]-[tau](P,Q))], then there exists a pair ([alpha], [beta]) of a non-negative number a and a real number [beta]such that

(2.1) [[sigma].sub.min] [P - (a + [[beta].sub.i])Q] = [[sigma].sub.min] [P - ([alpha] + [[eta].sub.i] + [[beta].sub.i])Q] = [delta].

We omit the proof of this theorem, since it is similar to that of [8, Theorem 3.1] but with horizontal search replaced by vertical search. If we rewrite problem (1.5) as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

the horizontal search can also be applied to solve (1.5). For the numerical verification, we only need a relation implied by (2.1).

COROLLARY 2.2. Assume that [delta] > [tau](P,Q). Given a number [eta] [member of] (0, 2([delta]- [tau](|P,Q))/[parallel]Q[[parallel].sub.2]], then there exists a pair ([alpha], [beta]) of a non-negative number a and a real number [beta] such that

(2.2) [delta] [member of] [sigma] [P - ([alpha] + [[beta].sub.i])Q] [intersection] [sigma] [P - ([alpha] + [[eta].sub.i] + [[beta].sub.i])Q],

where [alpha](*) is the set of all singular values of its matrix argument.

For [[delta].sub.1] > [[delta].sub.2] > 0, set [delta] = [[delta].sub.1] and [eta] = 2([[delta].sub.1] - [[delta].sub.2]/[parallel]Q[[parallel].sub.2]). This corollary implies that when no pair ([alpha], [beta]) satisfying (2.2) exists, the inequality [eta] > 2([delta]- [tau](P,Q))/[parallel]Q[[parallel].sub.2] is valid, so condition

(2.3) [tau](P,Q) > [[delta].sup.2]

holds. On the other hand, when a pair exists, then by definition we can conclude that

(2.4) [eta](P,Q) [less than or equal to] [[delta].sub.1].

This is called the modified Gu test for the stabilizability radius.

3. Trisection algorithm for computing stabilizability radii. The bisection algorithm of Gu  keeps only an upper bound on the distance to uncontrollability. It refines the upper bound until condition (2.3) is satisfied, and at termination the controllability radius lies within a factor of 2 of [[delta].sub.1], i.e., [[delta].sub.1]/2 [tau] t(P, Q) [less than or equal to] 2[[delta].sub.1]. To obtain the distance to uncontrollability with better accuracy, Burke et al.  proposed a trisection variant. Since the derivation of the details of the algorithm presented in this section follows  step by step, we just give brief results. Using the modified Gu test for the stabilizability radius instead of Gu's test, the trisection algorithm yields bounds of t(P, Q) in form of an interval [l, u] and reduces the length of this interval by a factor of 2/3 at each iteration.
```Trisection algorithm for computing (1.5)

Input:    P, Q [member of] [C.sup.kxn] with k [less than or equal to]
n and rank(Q) = k and a tolerance [epsilon] > 0.

Output:   Scalars l and u satisfying l < [tau](P, Q) [less than or
equal to] u and u - l < [epsilon].

Initialize the lower bound as l = 0 and the upper bound as u =
[[sigma].sub.min] (P).
repeat
[[delta].sub.1] = l + 2/3 (u - l)

[[delta].sub.2] = l + 1/3 (u - l)

Apply the modified Gu test for stabilizability radius.
if (2.4) is verified then
u [left arrow] [[delta].sub.1].
else
% Otherwise (2.3) is verified.
l [left arrow] 52.
end if
until u - l < [epsilon]
Return l and u.
```

3.1. Verification scheme. Let the singular value decomposition of Q be

Q = U [[SIGMA] 0] V*,

and define

[[P.sub.1] [P.sub.2]] := U*PV,

where [P.sub.1] [member of] [C.sup.pxp] and [P.sup.2] [member of] [C.sup.pxn-p]. It is remarkable that [delta] [member of] [sigma][P - ([alpha] + [[beta].sub.i])Q] if and only if a is an eigenvalue of the following matrix

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [[??].sub.2] := [P.sub.2][P.sup.*.sub.2]/[delta] - [delta]I. Then, condition (2.2) is equivalent to the fact that the following Sylvester equation

H ([beta],[delta])X + XH ([beta] + [eta],[delta]) =0

has a nontrivial solution. By the Kroneckerization of this Sylvester equation, for the verification of a pair ([alpha], [beta]) satisfying (2.2), we search at first for an imaginary eigenvalue i[beta] of the following matrix

(3.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where

[A.sub.1] := I [cross product] ([[SIGMA].sup.-1][P.sub.1]), [A.sub.2] := [([[SIGMA].sup.-1][P.sub.1] - i[eta]I).sup.T] [cross product] I,

[A.sub.3] := I [cross product] ([[SIGMA].sup.-1][P.sup.*.sub.1]), [A.sub.4] := [([[SIGMA].sup.-1][P.sup.*.sub.1] + i[eta]I).sup.T] [cross product] I,

[B.sub.1] := I [cross product] ([[SIGMA].sup.-1][[??].sub.2]), [B.sub.2] := [([[SIGMA].sup.-1][[??].sub.2]).sup.T] [cross product] I,

[C.sub.1] := I [cross product] ([[SIGMA].sup.-1]), [C.sub.2] := [([[SIGMA].sup.-1]).sup.T] [cross product] I.

Next, if there exists an imaginary eigenvalue i[beta] of (3.1) such that the matrices H([beta], [delta]) and H([beta] + [eta], [delta]) share a common non-negative eigenvalue [alpha], then the verification succeeds.

3.2. Eigenvalue search. We can observe that searching for an imaginary eigenvalue i[beta] of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the key operation of the modified Gu test for the stabilizability radius. In , the Matlab function eig can solve the search problem at a cost of O([k.sup.6]). In , an algorithm for searching real eigenvalues of a matrix in [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is given at a cost of O([p.sup.4]) on average by employing a shift-and-invert Arnoldi method of the matrix at a cost of O([p.sup.3]). Clearly, we can search for a real eigenvalue [beta] of - iA instead of searching for an imaginary eigenvalue i[beta] of A. Moreover, the linear system (- iA - vI) v = u can be solved at a cost of O([p.sup.3]) for any v [member of] R, by means of a Sylvester equation solver (such as the LAPACK routine dtrsyl ) as indicated by the following lemma.

LEMMA 3.1. Let u [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and suppose that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the solution of the linear system (-iA - vI)v = iu. Then the following Sylvester equation

(3.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

has a unique solution Z = [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] where u = [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and v = [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

4. Numerical experiments. Setting P := A and Q := i, we observe that [tau](P, Q) is the stability radius of the system (1.1). Hence, we can compare the accuracy of our method given in this paper and the method for computing the stability radius given in [2, 13]. Table 4.1 presents this comparison for a variety of examples chosen from EigTool  with a tolerance of 10-4. All the tests are run using Matlab 6.5 under Linux on a PC.

For the next results, we again choose matrices A from EigTool  and let B be matrices with normally distributed entries. Table 4.2 displays the results of [[tau].sub.0](A, B) and [[tau].sub.1](A) in the second and the third columns, respectively, when we apply our method for computing problem (1.5). These results illustrate the inequality given in  that to (A, B) [less than or equal to] [[tau].sub.1](A). In particular, the result for the Godunov and Skew Laplacian matrices show the importance of computing the stabilizability radius when only the system matrix A is perturbed.

Acknowledgment. The Matlab software used for computing the general problem (1.5) in the numerical experiments is based on software provided by Emre Mengi. The authors thank the referees for their careful reading of the manuscript and insightful comments.

REFERENCES

 E. ANDERSON, Z. BAI, C. BISCHOF, S. BLACKFORD, J. DEMMEL, J. DONGARRA, J. DU CROZ, A. GREENBAUM, S. HAMMARLING, A. MCKENNEY, AND D. SORENSON, LAPACK Users s Guide, 3rd ed., SIAM, Philadelphia, 1999.

 S. BOYD, V. BALAKRISHNAN, A regularity result for the singular values of a transfer matrix and a quadratically convergent algorithm for computing its Leo-norm, Systems Control Lett., 15 (1990), pp. 1-7.

 J. V. BURKE, A. S. LEWIS, AND M. L. OVERTON, Pseudospectral components and the distance to uncontrollability, SIAM J. Matrix Anal. Appl., 26 (2004), pp. 350-361.

 R. BYERS, Detecting nearly uncontrollable pairs, in Proceedings of the International Symposium MTNS89, Vol. III, M. A. Kaashoek, J. H. van Schuppen, and A. C. M. Ran, eds., Birkhauser, Boston, 1990, pp. 447-457.

 R. EISING, The distance between a system and the set of uncontrollable systems, in Mathematical Theory of Networks and Systems (Beer Sheva, 1983), P. A. Fuhrmann, ed., Lecture Notes in Control and Information Sciences, 58, Springer, London, 1984, pp. 303-314.

 --, Between controllable and uncontrollable, Systems Control Lett., 4 (1984), pp. 263-264.

 M. GAO AND M. NEUMANN, A global minimum search algorithm for estimating the distance to uncontrollability, Linear Algebra Appl., 188/189 (1993), pp. 305-350.

 M. GU, New methods for estimating the distance to uncontrollability, SIAM J. Matrix Anal. Appl., 21 (2000), pp. 989-1003.

 M. GU, E. MENGI, M. L. OVERTON, J. XIA, AND J. ZHU, Fast methods for estimating the distance to uncontrollability, SIAM J. Matrix Anal. Appl., 28 (2006), pp. 477-502.

 M. L. J. HAUTUS, Controllability and observability conditions of linear autonomous systems, Indag. Math., 31 (1969), pp. 443-448.

 C. HE, Estimating the distance to uncontrollability: a fast method and a slow one, Systems Control Lett., 26 (1995), pp. 275-281.

 D. C. KHANH AND D. D. X. THANH, Controllability radii and stabilizability radii of time-invariant linear systems, Vietnam J. Math., 34 (2006), pp. 495-499.

 E. MENGI, Measures for Robust Stability and Controllability, PhD Thesis, Department of Computer Science, Courant Institute of Mathematical Sciences, New York University, New York, 2006.

 C. C. PAIGE, Properties of numerical algorithms related to controllability, IEEE Trans. Automat. Control, 26 (1981), pp. 130-138.

 T. G. WRIGHT, EIGTOOL: A Graphical Tool for Nonsymmetric Eigenproblems, Oxford University Computing Laboratory, Oxford, 2002. http://www.comlab.ox.ac.uk/pseudospectra/eigtool/

D. C. KHANH ([dagger]), H. T. QUYEN ([double dagger]) and D. D. X. THANH ([section])

* Received January 13, 2009. Accepted August 30, 2013. Published online on December 9, 2013. Recommended by P. Van Dooren. This work was supported by the Vietnam National Foundation for Science and Technology Development Grant No. 101.01-2012.12.

([dagger]) HCM University of Technology, Dist. 10, Ho Chi Minh City, Vietnam (dckhanh@hcmhutech.edu. vn).

([double dagger]) HCM University of Architecture, Dist. 3, Ho Chi Minh City, Vietnam (thucquyen911@gmail .com).

([section]) John von Neumann Institute, Vietnam National University HCM, Thu Duc District, Ho Chi Minh City, Vietnam (thanh.duong@jvn.edu.vn).
```Table 4.1

Stability radii with tolerance = [10.sup.-4].

Example                        New method       Method in 

Airy(5)                    (0.00370, 0.00380]      0.0038
Airy(10)                   (0.01245, 0.01254]      0.0125
Convection Diffusion(5)    (0.60395, 0.60403]      0.6040
Convection Diffusion(10)   (0.75310, 0.75317]      0.7532
Transient(5)               (0.02935, 0.02942]      0.0294
Transient(10)              (0.02025, 0.02032]      0.0203

Table 4 2
Stabilizability radii with tolerance = [10.sup.-4].

Example                      [[tau].sub.0] (A,B)

Airy(5,2)                     (0.03777, 0.03786]
Airy(10,4)                    (0.16442, 0.16449]
Basor Morrison(5,2)           (0.72733, 0.72740]
Basor Morrison(10,4)          (0.76689, 0.76696]
Chebyshev(5,2)                (0.75035, 0.75044]
Chebyshev(10,4)               (0.82703, 0.82711]
Companion(5,2)                (0.42431, 0.42438]
Companion(10,4)               (0.65862, 0.65868]
Convection Diffusion(5,2)     (0.76432, 0.76439]
Convection Diffusion(10,4)    (2.08213, 2.08222]
Davies(5,2)                   (0.23170, 0.23176]
Davies(10,4)                  (0.70003, 0.70012]
Demmel(5,2)                   (0.10152, 0.10159]
Demmel(10,4)                  (0.27415, 0.27424]
Frank(5,2)                    (0.45907, 0.45916]
Frank(10,4)                   (0.69384, 0.69393]
Gallery(5,2)                  (0.50544, 0.50551]
Gauss Seidel(5,2)             (0.06383, 0.06392]
Gauss Seidel(10,4)            (0.07038, 0.07046]
Godunov(7,3)                  (1.39651, 1.39659]
Grcar(5,2)                    (0.49571, 0.49579]
Grcar(10,4)                   (0.44178, 0.44185]
Hatano(5,2)                   (1.31802, 1.31809]
Hatano(10,4)                  (0.44865, 0.44872]
Kahan(5,2)                    (0.18594, 0.18601]
Kahan(10,4)                   (0.05592, 0.05599]
Landau(5,2)                   (0.41766, 0.41773]
Landau(10,4)                  (0.36615, 0.36623]
Markov Chain(5,2)             (0.04355, 0.04365]
Markov Chain(10,4)            (0.20548, 0.20557]
Orr Sommerfield(5,2)          (0.04834, 0.04841]
Orr Sommerfield(10,4)         (0.08401, 0.08408]
Skew Laplacian(8,3)          (27.86320, 27.86330]
Supg(4,2)                     (0.06546, 0.06554]
Supg(9,4)                     (0.03810, 0.03817]
Transient(5,2)                (0.38264, 0.38273]
Transient(10,4)               (0.34967, 0.34974]
Twisted(5,2)                  (0.14929, 0.14936]
Twisted(10,4)                 (0.77197, 0.77204]

Example                       [[tau].sub.1] (A)

Airy(5,2)                     (0.04286, 0.04296]
Airy(10,4)                    (0.16822, 0.16829]
Basor Morrison(5,2)           (1.45174, 1.45181]
Basor Morrison(10,4)          (1.41940, 1.41948]
Chebyshev(5,2)                (1.61063, 1.61073]
Chebyshev(10,4)               (2.50581, 2.50591]
Companion(5,2)                (0.71518, 0.71525]
Companion(10,4)               (0.75474, 0.75481]
Convection Diffusion(5,2)     (1.12964, 1.12971]
Convection Diffusion(10,4)    (3.10039, 3.10048]
Davies(5,2)                   (5.51655, 5.51663]
Davies(10,4)                  (6.19431, 6.19440]
Demmel(5,2)                   (0.10152, 0.10159]
Demmel(10,4)                  (0.27417, 0.27425]
Frank(5,2)                    (0.80848, 0.80858]
Frank(10,4)                   (0.85247, 0.85255]
Gallery(5,2)                  (0.43033, 0.43042]
Gauss Seidel(5,2)             (0.06323, 0.06332]
Gauss Seidel(10,4)            (0.07039, 0.07046]
Godunov(7,3)                 (425.81455, 425.81464]
Grcar(5,2)                    (0.57320, 0.57328]
Grcar(10,4)                   (0.51903, 0.51911]
Hatano(5,2)                   (0.48167, 0.48175]
Hatano(10,4)                  (0.67699, 0.67708]
Kahan(5,2)                    (0.18710, 0.18717]
Kahan(10,4)                   (0.05651, 0.05658]
Landau(5,2)                   (1.04564, 1.04572]
Landau(10,4)                  (0.34016, 0.34023]
Markov Chain(5,2)             (0.04470, 0.04477]
Markov Chain(10,4)            (0.20852, 0.20861]
Orr Sommerfield(5,2)          (0.05296, 0.05305]
Orr Sommerfield(10,4)         (0.09904, 0.09911]
Skew Laplacian(8,3)          (32.94072, 32.94079]
Supg(4,2)                     (0.06612, 0.06619]
Supg(9,4)                     (0.03638, 0.03644]
Transient(5,2)                (0.39787, 0.39796]
Transient(10,4)               (0.35715, 0.35722]
Twisted(5,2)                  (0.20265, 0.20273]
Twisted(10,4)                 (1.01131, 1.01140]
```
COPYRIGHT 2013 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.