Preliminary Results on Generalized Fixed Order Interpolation - - PowerPoint PPT Presentation

Preliminary Results on Generalized Fixed Order Interpolation

Constantino Lagoa1

1Electrical Engineering Department

The Pennsylvania State University USA

Workshop on Uncertain Dynamical Systems 2011 Work done in collaboration with Chao Feng and Mario Sznaier

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 1 / 27

A General Interpolation Problem

Given Set of pairs of complex numbers F . = {(z0, F0), (z1, F1), . . . , (zNNP, FNNP)}, Set of real values T . = {g0, g1, . . . , gNCF} Supply function s(u, y) Bound on system order N.

Interpolation Problem

Find a causal SISO system G(z) of order no greater N satisfying G(zk) = Fk, k = 0, . . . , NNP; G(z) = g0 + g1z−1 + . . . + gNCFz−NCF + . . .; G(z) is passive with given supply function s(u, y).

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 2 / 27

A General Interpolation Problem

Given Set of pairs of complex numbers F . = {(z0, F0), (z1, F1), . . . , (zNNP, FNNP)}, Set of real values T . = {g0, g1, . . . , gNCF} Supply function s(u, y) Bound on system order N.

Interpolation Problem

Find a causal SISO system G(z) of order no greater N satisfying G(zk) = Fk, k = 0, . . . , NNP; G(z) = g0 + g1z−1 + . . . + gNCFz−NCF + . . .; G(z) is passive with given supply function s(u, y).

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 2 / 27

A General Interpolation Problem

Given Set of pairs of complex numbers F . = {(z0, F0), (z1, F1), . . . , (zNNP, FNNP)}, Set of real values T . = {g0, g1, . . . , gNCF} Supply function s(u, y) Bound on system order N.

Interpolation Problem

Find a causal SISO system G(z) of order no greater N satisfying G(zk) = Fk, k = 0, . . . , NNP; G(z) = g0 + g1z−1 + . . . + gNCFz−NCF + . . .; G(z) is passive with given supply function s(u, y).

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 2 / 27

A General Interpolation Problem

Given Set of pairs of complex numbers F . = {(z0, F0), (z1, F1), . . . , (zNNP, FNNP)}, Set of real values T . = {g0, g1, . . . , gNCF} Supply function s(u, y) Bound on system order N.

Interpolation Problem

Find a causal SISO system G(z) of order no greater N satisfying G(zk) = Fk, k = 0, . . . , NNP; G(z) = g0 + g1z−1 + . . . + gNCFz−NCF + . . .; G(z) is passive with given supply function s(u, y).

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 2 / 27

Application to Control Systems

This is a general framework that includes many problems in the control area Fixed order H∞ controller design. Simultaneous stabilization. Fixed order system system identification. Spectral estimation.

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 3 / 27

Related Work

Nevalinna-Pick interpolation Caratheodory-Fejer interpolation Mixed-domain system identification Polynomial approaches to fixed-order system identification Many more.....

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 4 / 27

Our Contribution

The general interpolation problem formulated here is a very complex problem.

Objective

Develop computationally efficient relaxations of this problem by Reformulating the problem as finding a point in a semialgebraic set Using results from the areas of polynomial optimization and sparsification

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 5 / 27

Our Contribution

The general interpolation problem formulated here is a very complex problem.

Objective

Develop computationally efficient relaxations of this problem by Reformulating the problem as finding a point in a semialgebraic set Using results from the areas of polynomial optimization and sparsification

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 5 / 27

Setup

Fixed order plant G(z) = b0 + b1z−1 + · · · + bmz−m 1 + a1z−1 + · · · + anz−n Design parameters Θ . = [−a1, . . . , −an, b0, . . . , bm]T Supply function s(u, y) . =

• u

y π11 π12 π12 π22 u y

• Lagoa (Penn State)

Fixed Order Interpolation Udine 2011 6 / 27

Consistency Set as a Semialgebraic Set

Time Domain Interpolation

Time domain interpolation conditions are linear equalities on the decision variables.

Algebraic Description

The time domain interpolation conditions are satisfied if and only if pCF

k (Θ) .

= ΘT · φk − gk = 0, k = 0, . . . , NCF where φk . = [gk−1, . . . , gk−n, uk, . . . , uk−m]T is the regressor and uk is the impulse function.

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 7 / 27

Consistency Set as a Semialgebraic Set

Frequency Domain Interpolation

Frequency domain conditions are also linear equalities on the decision variables.

Algebraic Description

The frequency domain interpolation conditions are satisfied if and only if pNP

2k (Θ)

. = Re{b(zk) − Fka(zk)} = 0, pNP

2k+1(Θ)

. = Im{b(zk) − Fka(zk)} = 0 for k = 0, . . . , NNP.

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 8 / 27

Consistency Set as a Semialgebraic Set

Stability - Jury’s Criterion

Jury’s array an an−1 an−2 · · · a2 a1 1 1 a1 a2 · · · an−2 an−1 an J1,n−1 J1,n−2 J1,n−3 · · · J1,1 J1,0 J1,0 J1,1 J1,2 · · · J1,n−2 J1,n−1 J2,n−2 J2,n−3 J2,n−4 · · · J2,0 J2,0 J2,1 J2,2 · · · J2,n−2 . . . . . . . . . Jn−2,2 Jn−2,1 Jn−2,0 Jk,i . =

• Jk−1,n−k

Jk−1,i 1 Jk−1,n−k−i

• ,

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 9 / 27

Consistency Set as a Semialgebraic Set

Stability - Jury’s Criterion

Lemma

The system G(z) is stable, i.e., all the roots of a(z) locate inside the unit circle, if and

• nly if the following inequalities hold,

1 +

n

• i=1

ai > 1 +

n

• i=1

(−1)iai > |an| < 1 |Jk,n−k| < |Jk,0|, 1 ≤ k ≤ n − 2 where Jk,i, k = 1, . . . , n − 2, i = 0, . . . , n − k are the elements in the Jury’s array.

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 10 / 27

Consistency Set as a Semialgebraic Set

Stability - Jury’s Criterion

Hence, stability is equivalent to satisfaction of polynomial inequalities.

Algebraic Description

The system G(z) is stable if and only if pS

k (Θ) > 0, k = 1, . . . , n + 2.

where pS

1 (Θ)

. = 1 +

n

• i=1

ai, pS

2 (Θ)

. = 1 +

n

• i=1

(−1)iai, pS

3 (Θ)

. = 1 + an, pS

4 (Θ)

. = 1 − an, pS

k+4(Θ)

. = J2

k,0 − J2 k,n−k, k = 1, . . . , n − 2.

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 11 / 27

Consistency Set as a Semialgebraic Set

Dissipativity - KYP Lemma

Theorem

Assume that G(z) is stable and has no poles on the unit circle, then it is dissipative with respect to the quadratic supply function s(ui, yi) . = ui yi T Π ui yi

• ,

where Π ∈ R2×2 is a symmetric matrix, if and only if

• 1

G(z) ∗ Π

• 1

G(z)

• ≤ 0

for all z ∈ C, |z| = 1.

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 12 / 27

Consistency Set as a Semialgebraic Set

Dissipativity - Positive Polynomials on the Unit Circle

Note that

• 1

G(z) ∗ Π

• 1

G(z)

• ≤ 0

⇔ R(z) ≥ 0 where R(z) is the trigonometric polynomial R(z) . = −

• π11a(z)a(z−1) + π12a(z)b(z−1) + π12b(z)a(z−1) + π22b(z)b(z−1)
• Lemma

A trigonometric polynomial R(z) is non-negative on the unit circle if and only if there exists a dr + 1 × dr + 1 positive semi-definite matrix Q . = [qi,j] 0. such that rk =

dr +k

• i=k

qi,i−k, k = 0, . . . , dr.

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 13 / 27

Consistency Set as a Semialgebraic Set

Dissipativity - Positive Polynomials on the Unit Circle

Note that

• 1

G(z) ∗ Π

• 1

G(z)

• ≤ 0

⇔ R(z) ≥ 0 where R(z) is the trigonometric polynomial R(z) . =

dr

• k=−dr

rizi. where

✄ ✂

• ✁

ri = r−i

Lemma

A trigonometric polynomial R(z) is non-negative on the unit circle if and only if there exists a dr + 1 × dr + 1 positive semi-definite matrix Q . = [qi,j] 0. such that rk =

dr +k

• i=k

qi,i−k, k = 0, . . . , dr.

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 13 / 27

Consistency Set as a Semialgebraic Set

Dissipativity - Positive Polynomials on the Unit Circle

Note that

• 1

G(z) ∗ Π

• 1

G(z)

• ≤ 0

⇔ R(z) ≥ 0 where R(z) is the trigonometric polynomial R(z) . =

dr

• k=−dr

rizi. where

✄ ✂

• ✁

ri = r−i

Lemma

A trigonometric polynomial R(z) is non-negative on the unit circle if and only if there exists a dr + 1 × dr + 1 positive semi-definite matrix Q . = [qi,j] 0. such that rk =

dr +k

• i=k

qi,i−k, k = 0, . . . , dr.

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 13 / 27

Consistency Set as a Semialgebraic Set

Dissipativity - Algebraic formulation

Parametrization of all Q with additional free variables. For n = 2 and m = 1 ⇒ Q =   λ1 λ3 r2 ∗ λ2 r1 − λ1 ∗ ∗ r0 − λ1 − λ2   . Q is polynomial in Θ and λ

Lemma

Given a stable G(z) with no poles on the unit circle, G(z) is dissipative with respect to the quadratic supply matrix Π if and only if pDP

k (Θ, λ) ≥ 0, k = 1, . . . , dr + 1

for some λ where pDP

k (Θ, λ) are the principal minors of Q.

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 14 / 27

Consistency Set as a Semialgebraic Set

Dissipativity - Algebraic formulation

Parametrization of all Q with additional free variables. For n = 2 and m = 1 ⇒ Q =   λ1 λ3 r2 ∗ λ2 r1 − λ1 ∗ ∗ r0 − λ1 − λ2   . Q is polynomial in Θ and λ

Lemma

Given a stable G(z) with no poles on the unit circle, G(z) is dissipative with respect to the quadratic supply matrix Π if and only if pDP

k (Θ, λ) ≥ 0, k = 1, . . . , dr + 1

for some λ where pDP

k (Θ, λ) are the principal minors of Q.

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 14 / 27

Consistency Set as a Semialgebraic Set

Dissipativity - Algebraic formulation

Parametrization of all Q with additional free variables. For n = 2 and m = 1 ⇒ Q =   λ1 λ3 r2 ∗ λ2 r1 − λ1 ∗ ∗ r0 − λ1 − λ2   . Q is polynomial in Θ and λ

Lemma

Given a stable G(z) with no poles on the unit circle, G(z) is dissipative with respect to the quadratic supply matrix Π if and only if pDP

k (Θ, λ) ≥ 0, k = 1, . . . , dr + 1

for some λ where pDP

k (Θ, λ) are the principal minors of Q.

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 14 / 27

Consistency Set as a Semialgebraic Set

Summary

The consistency set is K . =                                    Θ : pCF

k (Θ) = 0,

k = 0, . . . , NCF pNP

k (Θ) = 0,

k = 0, . . . , 2NNP + 1 pS

k (Θ) ≥ 0,

k = 1, . . . , 2n pDP

k (Θ, λ) ≥ 0,

k = 1, . . . , dr + 1 for some λ                                   

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 15 / 27

SDP Relaxation

Find truncated moment sequence m satisfying MN(m) 0 MNi (pCF

k m) = 0,

k = 0, . . . , NCF MNi (pNP

k m) = 0,

k = 0, . . . , 2NNP + 1 MNi (pS

k m) 0,

k = 1, . . . , n + 2 MNi (pDP

k m) 0,

k = 1, . . . , dr + 1

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 16 / 27

SDP Relaxation

Moment Matrices

Example of moment matrix MN(m) =      M0,0(m) M0,1(m) · · · M0,N(m) M1,0(m) M1,1(m) · · · M1,N(m) . . . . . . ... . . . MN,0(m) MN,1(m) · · · MN,N(m)      Mj,k(m) =      mj+k,0 mj+k−1,1 · · · mj,k mj+k−1,1 mj+k−2,2 · · · mj−1,k+1 . . . . . . ... . . . mk,j mk−1,j+1 · · · m0,j+k      Localizing matrix MNi (pi, m) is defined as MNi (pi, m)(i, j) =

• α

pi,αm(β(i, j) + α)

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 17 / 27

SDP Relaxation

Moment Matrices

Example of moment matrix MN(m) =      M0,0(m) M0,1(m) · · · M0,N(m) M1,0(m) M1,1(m) · · · M1,N(m) . . . . . . ... . . . MN,0(m) MN,1(m) · · · MN,N(m)      Mj,k(m) =      mj+k,0 mj+k−1,1 · · · mj,k mj+k−1,1 mj+k−2,2 · · · mj−1,k+1 . . . . . . ... . . . mk,j mk−1,j+1 · · · m0,j+k      Localizing matrix MNi (pi, m) is defined as MNi (pi, m)(i, j) =

• α

pi,αm(β(i, j) + α)

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 17 / 27

SDP Relaxation

Relaxation converges to solution as N → ∞. However, complexity grows fast with N. We know that we have a solution if moment matrix is sparse i.e., low rank MN

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 18 / 27

SDP Relaxation

Relaxation converges to solution as N → ∞. However, complexity grows fast with N. We know that we have a solution if moment matrix is sparse i.e., low rank MN

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 18 / 27

SDP Relaxation

Relaxation converges to solution as N → ∞. However, complexity grows fast with N. We know that we have a solution if moment matrix is sparse i.e., low rank MN

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 18 / 27

Alternative Algorithm

m∗ = arg min

m

rank(MN(m)) s.t. MN(m) 0,

• α

pCF

k,αmα = 0, k = 0, . . . , NCF,

• α

pNP

k,αmα = 0, k = 0, . . . , 2NNP + 1,

• α

pS

k,αmα ≥ 0, k = 1, . . . , n + 2,

• α

pDP

k,αmα ≥ 0, k = 1, . . . , dr + 1,

N = dmax/2

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 19 / 27

Alternative Algorithm

m∗ = arg min

m

rank(MN(m)) s.t. MN(m) 0,

• α

pCF

k,αmα = 0, k = 0, . . . , NCF,

• α

pNP

k,αmα = 0, k = 0, . . . , 2NNP + 1,

• α

pS

k,αmα ≥ 0, k = 1, . . . , n + 2,

• α

pDP

k,αmα ≥ 0, k = 1, . . . , dr + 1,

✞ ✝ ☎ ✆

N = dmax/2

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 19 / 27

“Convexification” of Rank

Fazel et al.

Set X ← M(m), X0 ← I, δ ← 0, k ← 0. repeat Solve Xk+1 ← arg min Tr(Xk + δI)−1X s.t. LMI constraints. Decompose the symmetric matrix Xk = T −1DT. Set δ ← min diag(D) + δ0. Set k ← k + 1. until a convergence criterion is reached. return Xk

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 20 / 27

Remarks on System Identification

Noisy time measurements yk = gk + et

k;

|et

k| ≤ ¯

et Noisy frequency measurements Fk = G(e−jωk ) + ef

k;

|ef

k| ≤ ¯

ef

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 21 / 27

Remarks on System Identification

New frequency domain constraints |G(e−jωk ) − Fk| ≤ ¯ ef,

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 22 / 27

Remarks on System Identification

New time domain constraints ΘTφk + ΘT∆ηk = yk − et

k,

¯ et + et

k ≥ 0,

¯ et − et

k ≥ 0,

∆ηk: vector that depends (linearly) on the noise variables et

k−1, . . . , et k−n

Need new variables et

k and no longer linear constraints

Complexity grows linearly with number of time measurements Consequence of sparsity of constraints (Running intersection property)

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 23 / 27

Remarks on System Identification

New time domain constraints ΘTφk + ΘT∆ηk = yk − et

k,

¯ et + et

k ≥ 0,

¯ et − et

k ≥ 0,

∆ηk: vector that depends (linearly) on the noise variables et

k−1, . . . , et k−n

Need new variables et

k and no longer linear constraints

Complexity grows linearly with number of time measurements Consequence of sparsity of constraints (Running intersection property)

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 23 / 27

Remarks on System Identification

New time domain constraints ΘTφk + ΘT∆ηk = yk − et

k,

¯ et + et

k ≥ 0,

¯ et − et

k ≥ 0,

∆ηk: vector that depends (linearly) on the noise variables et

k−1, . . . , et k−n

Need new variables et

k and no longer linear constraints

Complexity grows linearly with number of time measurements Consequence of sparsity of constraints (Running intersection property)

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 23 / 27

Remarks on Fixed-order Controller Design

Gcl(z) = P0K[1 + P0K]−1 Interpolation conditions: Gcl(z) stable Gcl(pk) = 1 at unstable poles pk of P0 Gcl(zk) = 0 at unstable zeros zk of P0. Gcl∞ ≤ γ

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 24 / 27

Numerical Example 1

Stable Fixed Order System Identification

G(z) = z2 − 3z + 2 z3 − 0.6300z2 − 0.7228z + 0.8043. System has G(z−1)∞ = 58.32 and root radius ρ = 0.97 Random input uniformly distributed in [−1, 1]. Magnitude of noise 0.2 Attempt 1 - Data given to algorithm G(z−1)∞ ≤ 50 and root radius ρ ≤ 0.965 Input and bound on noise.

✞ ✝ ☎ ✆

Algorithm concluded consistency set is empty

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 25 / 27

Numerical Example 1

Stable Fixed Order System Identification

G(z) = z2 − 3z + 2 z3 − 0.6300z2 − 0.7228z + 0.8043. System has G(z−1)∞ = 58.32 and root radius ρ = 0.97 Random input uniformly distributed in [−1, 1]. Magnitude of noise 0.2 Attempt 1 - Data given to algorithm G(z−1)∞ ≤ 50 and root radius ρ ≤ 0.965 Input and bound on noise.

✞ ✝ ☎ ✆

Algorithm concluded consistency set is empty

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 25 / 27

Numerical Example 1

Stable Fixed Order System Identification

G(z) = z2 − 3z + 2 z3 − 0.6300z2 − 0.7228z + 0.8043. System has G(z−1)∞ = 58.32 and root radius ρ = 0.97 Random input uniformly distributed in [−1, 1]. Magnitude of noise 0.2 Attempt 1 - Data given to algorithm G(z−1)∞ ≤ 50 and root radius ρ ≤ 0.965 Input and bound on noise.

✞ ✝ ☎ ✆

Algorithm concluded consistency set is empty

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 25 / 27

Numerical Example 1

Stable Fixed Order System Identification

G(z) = z2 − 3z + 2 z3 − 0.6300z2 − 0.7228z + 0.8043. System has G(z−1)∞ = 58.32 and root radius ρ = 0.97 Random input uniformly distributed in [−1, 1]. Magnitude of noise 0.2 Attempt 2 - Data given to algorithm G(z−1)∞ ≤ 60 and root radius ρ ≤ 0.97 Input and bound on noise. ˆ G(z) = 0.9955z2 − 2.9690z + 1.9690 z3 − 0.605z2 − 0.7277z + 0.771

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 25 / 27

Numerical Example 1

Stable Fixed Order System Identification

G(z) = z2 − 3z + 2 z3 − 0.6300z2 − 0.7228z + 0.8043. System has G(z−1)∞ = 58.32 and root radius ρ = 0.97 Random input uniformly distributed in [−1, 1]. Magnitude of noise 0.2 Attempt 2 - Data given to algorithm G(z−1)∞ ≤ 60 and root radius ρ ≤ 0.97 Input and bound on noise. ˆ G(z) = 0.9955z2 − 2.9690z + 1.9690 z3 − 0.605z2 − 0.7277z + 0.771

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 25 / 27

Numerical Example 2

Stable Controller Design

Plant: P(z) = z(z − 3) (z − 2)(z + 2) Controller: F and K(z) = b1 + b0z a1 + z Objectives: Closed loop system Gc(z) = GK[1 + GK]−1 is internally stable; DC gain from r to y is 1; Controller K(z) is stable; Transfer function from r to u has H∞ norm less than 2.

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 26 / 27

Numerical Example 2

Stable Controller Design

Plant: P(z) = z(z − 3) (z − 2)(z + 2) Controller: F and K(z) = b1 + b0z a1 + z Controller obtained: F = 0.0928 K(z) = −0.2361z + 1.417 z

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 26 / 27

Concluding Remarks

In this presentation: Framework for fixed order interpolation Allows for conditions on

Time response Frequency response Stability Dissipativity

Proposed algorithm based on results from the areas of polynomial optimization and sparsity. Further research Improving numerical performance of algorithms Extending framework to other types of constraints

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 27 / 27

Concluding Remarks

In this presentation: Framework for fixed order interpolation Allows for conditions on

Time response Frequency response Stability Dissipativity

Proposed algorithm based on results from the areas of polynomial optimization and sparsity. Further research Improving numerical performance of algorithms Extending framework to other types of constraints

Lagoa (Penn State) Fixed Order Interpolation Udine 2011 27 / 27