Level set methods for robustness measures P. V AN D OOREN CESAME, - - PowerPoint PPT Presentation

Level set methods for robustness measures

• P. VAN DOOREN

CESAME, Universit´ e catholique de Louvain joint work with Y. Genin, T. Lawrence, M. Overton, J. Sreedhar, A. Tits

What it is about

Fast and reliable computation of robustness measures for

• Stability
• Passivity
• Minimality

for (possibly structured) perturbations {∆A, ∆B, ∆C, ∆D}

• f a (real or complex) plant

{A, B, C, D}

Link to structured singular values

In each of these problems we need to find a point λ ∈ C such that some singular value or eigenvalue of matrices derived from

• A − λI
• A − λI

B

• 

 A − λI C     A − λI B C D   is minimal. Level sets (parameterized by ξ) are sets where singular values or eigenvalues are larger than the parameter ξ. This is clearly related to so-called pseudo-spectra

Setting

Given a minimal system G(λ) := C(λI − A)−1B + D and a perturbed system G∆(λ) := C∆(λI − A∆)−1B∆ + D∆ we measure perturbations using ∆ :=   ∆A ∆B ∆C ∆D   :=   A∆ B∆ C∆ D∆   −   A B C D   and use the ∆2 norm All matrices can be real or complex

Complex stability radius

We are looking for the radius rC = inf∆{∆2 | A∆ is unstable} (Hinrichsen-Pritchard) We need to find an eigenvalue λ of A∆ in the unstable region Γ or on its boundary ∂Γ: det(A + ∆A − λI) = 0, λ ∈ ∂Γ For the complex case, perturbation theory says ∆A2 = σmin(A − λI) so we have (cfr pseudospectrum) rC = min

λ∈∂Γ σmin(A − λI) = 0,

λ ∈ ∂Γ

Resolvant

In systems theory we are more familiar with the resolvent and its norm : G(λ) := (λI − A)−1, r−1

C

= max

λ∈∂Γ σmax[G(λ)]

−2 −1 1 2 −2 −1.5 −1 −0.5 0.5 1 1.5 2

Complex stability radius

We wanted the radius rC(A) = inf∆∈Cm×p{∆2 | A + ∆A is unstable} We needed to find an eigenvalue λ of A + ∆A in the unstable region Γ or

• n its boundary ∂Γ:

det(A + ∆A − λI) = 0, λ ∈ ∂Γ For the complex case, this depended on the transfer function G(λ) = (λI − A)−1 and yielded r−1

C

= max

λ∈∂Γ σmax[G(λ)]

Structured stability radius

If we impose a very simple structure, the problem becomes rC(A, B, C) := inf

∆∈Cm×p{∆2 : A + B∆C is stable}

Define G(λ) := C(λI − A)−1B then r−1

C (A, B, C) = maxλ∈∂Γ{σmax[G(λ)]}

See Hinrichsen Pritchard for more on this

Maximum frequency response

We thus need a reliable algorithm to find (λ ∈ ∂Γ = jω or ejω): max

ω

σmaxC(ejωI − A)−1B

−4 −3 −2 −1 1 2 3 4 4 6 8 10 12 14 16 Amplitude of frequency response H(ejω)

Alternative : gradient search but too many peaks may cost a lot and may be trouble- some The solution is to look at the ξ level set

Maximum frequency response

We thus need a reliable algorithm to find (λ ∈ ∂Γ = jω or ejω): max

ω

σmaxC(ejωI − A)−1B

−4 −3 −2 −1 1 2 3 4 4 6 8 10 12 14 16 Amplitude of frequency response H(ejω) ξ

Alternative : gradient search but too many peaks may cost a lot and may be trouble- some The solution is to look at the ξ level set

Bisection

The intersection points ωi with the ξ level (of all singular values) are imaginary eigenvalues of a Hamil- tonian (or symplectic) matrix H(ξ) :=   A −BB∗/ξ C∗C/ξ −A∗   This can be used to do e.g. bisection Linear convergence (Byers)

10 10

1

10

2

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Frequency ω2 ω4 ξ0 ≡ ξold ω6 ω8 ω7 ω3 ω5 ω1

Interval midpoint rule is quadratic (Boyd et al, Bruinsma et al)

Correct intervals and interpolation

The eigenvalue problem also yields the derivatives of the singular value plots They yield the info to find each indi- vidual singular value

−4 −3 −2 −1 1 2 3 4 2 4 6 8 10 12 Two dominant singular values frequency

Correct intervals and interpolation

The eigenvalue problem also yields the derivatives of the singular value plots They yield the info to find each indi- vidual singular value

−4 −3 −2 −1 1 2 3 4 2 4 6 8 10 12 Two dominant singular values frequency

Midpoint vs interpolation

Acceleration techniques yielding superlinear convergence

2 2.2 2.4 2.6 2.8 3 10 10.5 11 11.5 12 12.5 13 σmax midpoint rule fmid=12.76, fint=12.82 frequency

Midpoint rule (Boyd et al) Choose ω+ as midpoint of interval Choose ξ+ as function at that ω+ |ξ+ − ξ∗| = O|ξ − ξ∗|2 Interpolation rule (Genin-V) Choose ξ++ as maximum of cubic interpolating polynomial |ξ++ − ξ∗| = O|ξ − ξ∗|4

Quartic convergence

Iteration ξ (midp.) Intervals (midp.) ξ (cubic) Intervals (cubic) 1 0.5224 [0,1.1991] 0.5224 [0,1.1991] 2 0.7980 [0.1867,0.5995] [0.7097,1.0153] 6.5148 [0.7804,0.7994] 3 1.7669 [0.7472,0.8625] 8.4043 [0.78942,0.78943] 4 5.3027 [0.7762,0.8048] 8.4043 Convergence 5 8.3691 [0.7884,0.7905] 6 8.4043 [0.78942,0.78943]

Complexity results (per iteration): ωi : a(2n)3, a ≈ 50 (exploit Hamiltonian structure) ∂σj/∂ω : bn2, b < n (exploit Hamiltonian structure) ξk : cn2(m + p), c < n (use condensed forms) Overall complexity is O(n)3

Real stability radius

rR(A, B, C) := inf

∆∈Rm×p{∆2 : A + B∆C is stable}

Define G(λ) := C(λI − A)−1B and G := Gr + iGi then (Qiu et al) r−1

R (A, B, C) = supλ∈∂Γ{µR[G(λ)]}

where µR(G) := inf

γ∈(0,1] σ2 [Pγ] :=

inf

γ∈(0,1] σ2

  Gr −γGi Gi/γ Gr  

Associated Hamiltonian

We use Tγ :=

1 √ 2

  I γI I/γ −I   ,   Gr −γGi Gi/γ Gr   = Tγ   G G   Tγ to derive an associated real Hamiltonian matrix Hγ(ξ) :=        A αBBT −βBBT −A −βBBT αBBT −αCT C βCT C −AT βCT C −αCT C AT        where α := (1 + γ2)/(2γξ), β := (1 − γ2)/(2γξ) We can find the jω eigenvalue of Hγ(ξ) that correspond to σ2[Pγ(ω)]

Lower envelope

For each γo value there is a σ2 plot whose levels we can check with Hγo(ξo)

1 2 3 4 5 6 7 8 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 ω ξ λγ(H(ω)) λmax(Hγ

• (ω))

The (solid) curve µR(ω) we need to maximize is the lower envelope

• f these (dotted) curves

Each one is tangent to µR(ω) in

• ne frequency ωo

Convergence is quadratic or cubic (Sreedhar-Tits-V)

Passivity radius

Let G(λ) := C(λIn − A)−1B + D be strictly passive i.e. stable and positive real Re λi(A) < 0, G(jω) + [G(jω)]∗ ≻ 0, ∀ω ∈ R Consider the perturbed system G∆(λ) := C∆(λIn − A∆)−1B∆ + D∆ We wish to find the passivity radius of the system G(λ) prC(G) := inf

∆ {∆2 | G∆(λ) is not passive}.

KYP lemma

Passivity (stability and positive realness) Re λi(A) < 0, G∆(jω) + [G∆(jω)]∗ ≻ 0, ∀ω ∈ R iff there exists a Hermitian matrix P such that   −A∆P − PA∗

∆

B∆ − PC∗

∆

B∗

∆ − C∆P

D∆ + D∗

∆

  ≻ 0, P ≻ 0 Stability Re λi(A) < 0 iff there exists a Hermitian matrix P such that −A∆P − PA∗

∆ ≻ 0,

P ≻ 0 Therefore stability can not be lost “before” positive realness is lost (It can at one frequency if minimality is also lost)

Positive realness G∆(jω) + [G∆(jω)]∗ ≻ 0 gets lost as soon as det     A∆ − jωIn B∆ A∗

∆ + jωIn

C∗

∆

B∗

∆

C∆ D∆ + D∗

∆

    = 0, for some ω ∈ R

• r as soon as

det  Hω + E   ∆ ∆∗   ET   = 0 where Hω :=     A − jωIn B A∗ + jωIn C∗ B∗ C D + D∗     , E :=     In In Im Im    

We need a closed expression for min ∆2 : det  Hω + E   ∆ ∆∗   ET   = 0 The corresponding transfer function is H(ω) := ET H(ω)−1E and one shows (Hu-Qiu, Overton-V) that the passivity radius is then pr−1

C (A, B, C, D) = supω{νC[H(ω)]}

where νC := max{inf

γ λmax(Hγ), inf γ λmax(−Hγ)}

and Hγ := TγHTγ, Tγ :=   γIn+m In+m/γ  

Associated Hamiltonian

For each γo value there is a λmax plot whose levels we can check with   −γ2

• In/ξo

A − jωIn A∗ + jωIn −γ−2

• In/ξo

 −   B C∗   (D+D∗−γ2

• + γ−2
• ξo

I)−1 B∗ C

• 1

2 3 4 5 6 7 8 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 ω ξ λγ(H(ω)) λmax(Hγ

• (ω))

The (solid) curve is the lower enve- lope of the dotted curves and each

• ne is tangent to it in one frequency

ωo Convergence is quadratic or cubic (Overton-V)

Real passivity

An equivalent problem has been formulated (Hu-Qiu) but only a bound can be obtained in general It is still an open problem if this bound is always met for the passivity radius. This bound can be computed using a level set method involving two parameters γ1 and γ2 The optimization problem over γ1 and γ2 is simple (quasi-convex)

Minimality

A system remains controllable unless it is perturbed by the amount mrC(A, B) := min{[∆A | ∆B], (A+∆A, B+∆B) is uncontrollable } For complex perturbations this is non-convex minimization (Eising) mrC(A, B) := τC(A, B) := min

λ∈C σmin [A − λI | B]

Re(z) Im(z) σmin([A−zI,B]) −2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 2.5 −6 −5 −4 −3 −2 −1

Bounded derivative

The slope of the singular value plot is bounded : |∂σi ∂λ | = |u∗[In | 0]v| ≤ 1 We have the following (Gu) Theorem For any δ > τ(A, B) and for any η ∈ [0, 2(δ −τ)] there exist two pairs of real numbers α, β such that δ ∈ σ[A − (α + βi)I, B] δ ∈ σ[A − (α + η + βi)I, B]

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.2 0.4 0.6 0.8 1 Singular value plot and derivative bounds

τ

Bisection - trisection

The above theorem yields one of the following bounds (Mengi et al) δ ≥ τ,

• r

τ > (δ − η)/2 This yields a bisection (or trisection) algorithm to estimate τC(A, B) within a factor 2 at most. The algorithm is based on the solution of a Sylvester equation to test common roots of two eigenvalue problems The complexity is O(n6) for the basic algorithm but between O(n5) and O(n4) using sparse techniques No guaranteed accurate computation unless one can find isolated quasi-convex regions Can be extended to the estimation of τR(A, B)

Other uses of level sets

• An upper bound for µ (Lawrence-Tits-V)
• Radius of hyperbolicity and ellipticity of second order polynomial

matrices P(λ) := λ2A + λB + C (Hachez-V)

• Radius of definiteness of pencils λB − A (Crawford number)

(Higham-Tisseur-V) Flexible and fast method for problems with a characterization in terms of eigenvalues or singular values Depending on frequency and in a quasi-convex way of a few parameters

Some references with level set ideas

• R. BYERS, A bisection method for measuring the distance of a stable matrix to the unstable

matrices, SIAM J. Sci. Stat. Comput., vol. 9, pp. 875–881, Sept 1988.

• S. BOYD, V. BALAKRISHNAN, AND P. KABAMBA, A bisection method for computing the

H∞ norm of a transfer matrix and related problems, Mathematics of Control, Signals, and Systems, vol. 2, pp. 207–219, 1989.

• D. HINRICHSEN, B. KELB, AND A. LINNEMANN, An algorithm for the computation of the

complex stability radius, Automatica, vol. 25, pp. 771–775, 1989.

• N. BRUINSMA, AND M. STEINBUCH, A fast algorithm to compute the H∞ norm of a

transfer function matrix, Systems & Control Letters, vol. 14, pp. 287–293, 1990.

• J. SREEDHAR, P. VAN DOOREN, AND A.L. TITS, A fast algorithm to compute the real

structured stability radius, in Int. Series of Numer. Math., Birkhauser, pp. 219–230, 1996.

• Y. GENIN, P. VAN DOOREN, AND V. VERMAUT. Convergence of the calculation of H∞

norms and related questions. in Proceedings MTNS-98, pp. 429–432, July 1998.

• C. LAWRENCE, A. TITS, AND P. VAN DOOREN, A fast algorithm for the computation of an

upper bound on the µ-norm, Automatica, vol. 36, pp. 449–456, 2000.

• M. OVERTON, AND P. VAN DOOREN, On computing the complex passivity radius, in

Proceedings 44th IEEE Conference on Decision and Control, pp. 7960–7964, Spain, 2005.