ROBUST INPUT DESIGN USING SUM OF SQUARES CONSTRAINTS

14th IFAC Symposium on System Identification, Newcastle, Australia, 2006

ROBUST INPUT DESIGN USING SUM OF SQUARES CONSTRAINTS Jonas Mårtensson and Håkan Hjalmarsson Department of Signals, Sensors and Systems, KTH SE-100 44 Stockholm, Sweden [email protected]

Abstract: The objective of this paper is to illustrate how recent advances in convex programming can be used in input design in system identification. It is shown that sum of squares techniques and a recent generalization of the KYP lemma can be used to transform input design problems with frequency dependent bounds on the frequency function to convex optimization problems. Another important problem that can be handled with sum of squares is robust input design where initial uncertainty in the true parameters c is accounted for. Copyright 2006 IFAC Keywords: Input design, Sum of Squares, KYP-lemma.

1. INTRODUCTION

unknown system. The examples are solved with the toolbox YALMIP for MATLAB, see (Löfberg, 2004).

The input design problem has received significant interest recently (Hildebrand and Gevers, 2003), (Jansson and Hjalmarsson, 2005), (Bombois et al., 2004b) and (Bombois et al., 2004a). One driving force for this has been advances in convex programming. In this paper we continue this development and present some initial work where sum of squares (sos) constraints are used in the optimization problems and where we explore a recent generalization (Iwasaki and Hara, 2005) of the celebrated KYP-lemma (Yakubovich, 1962). A brief introduction to Sum of Squares is given in section 2. The input design problem we will use is formulated in section 3. The results are presented in three examples. In Section 4 we show how to recast frequency dependent, infinite dimensional, constraints over the entire frequency axis as finite dimensional sum of squares constraints. In Section 5, we show how to extend this to constraints defined on segments of the frequency axis, e.g. piecewise constant constraints. We also show that this type of problem can be handled more efficiently by the generalized KYP-lemma (Iwasaki and Hara, 2005). Optimal designs are system dependent and robustness with respect to the underlying system is an important issue (Hjalmarsson, 2004; Jansson and Hjalmarsson, 2005). A recent and interesting contribution to this area is (Goodwin et al., 2006). In Section 6 we show how sum of squares can be used to generate input designs that are robust with respect to the 1

This work was supported by the Swedish Research Council

1352

2. SUM OF SQUARES This section serves as a short introduction to the problem of finding global optima of multi-variable polynomial functions using sum of squares relaxations, see e.g. (Parrilo and Sturmfels, 2003). Consider a polynomial function f (x) of degree 2d in n variables, x ∈ Rn . The main idea is to exchange the positivity constraint f (x) ≥ 0 with ’ f (x) is a sum of squares’. Finding a sum of squares decomposition of a polynomial function is equivalent to solving a semi-definite program, (Parrilo, 2003). Express the polynomial f (x) as a quadratic form in all possible monomials of degree less than or equal to d, i.e. f (x) = zT Qz

(1)

where  Q is a constant matrix and the dimension of z is n+d d . If the matrix Q is positive semi-definite (psd) then there exist a decomposition Q = LT L and f (x) = ∑i (Lz)2i is a sum of squares. If f (x) has a sum of squares decomposition, then f (x) is positive. However, the representation (1) is not unique and Q may not be psd for some representations even though f (x) is positive. The problem of deciding whether a psd matrix Q that fulfills the constraints given by (1) exists, can be solved by the feasibility of a semidefinite program.

Λ G (θ ) Λ H (θ ) ∂ G(θ ) H(θ ) , Fe (θ ) = H(θ ) , ΛG (θ ) = ∂ θ θ) jω in = ∂ H( ∂ θ . The frequency arguments e

In the following sections, it will be shown how sum of squares constraints can be used to solve various input design problems, but first we need some introduction to the field of input design.

Fu (θ ) =

ΛH (θ ) e.g. G(e jω , θ ) are omitted for brevity. The complex conjugate transpose is denoted as (·)∗ . The expressions (3) and (4) show the influence on the parameter estimation error of different choices of input spectra and are therefore useful when considering input design. In open loop, the only quantity that can be used to affect the uncertainty in the parameters is actually the input spectrum Φu (ω ).

3. THE INPUT DESIGN PROBLEM The problem of interest is that of finding an optimal input signal for a system identification experiment. Of course, there are many different optimality conditions and many different signal constraints and quality constraints that can be used. In several contributions, the objective has been to optimize some performance measure subject to constraints on the input and/or output signals, see e.g. (Yuan and Ljung, 1985), (Gevers and Ljung, 1986), (Hjalmarsson et al., 1996) and (Zhu and van den Bosch, 2000). One could also consider optimizing some signal quantity subject to model quality constraints. Such approaches are taken in e.g. (Bombois et al., 2004b), (Jansson, 2004) and (Mårtensson et al., 2005).

Here, the main interest is not the accuracy of the parameter estimate θˆN , but rather the accuracy of the estimated transfer function G(e jω , θˆN ). Using a first order Taylor approximation of G around the true parameter value, the parameter uncertainty can be mapped to the uncertainty in the transfer function estimate as ˆ jω ) = ΛG (e jω , θo )∗ P(Φu , θo )ΛG (e jω , θo ), Var G(e (5) where the dependence of the input spectrum Φu and the true system θo is clearly visible in the notation.

In this paper we follow some of the ideas from the latter approach, but the main focus is not to emphasize a particular input design strategy. The purpose is to illustrate that sum of squares constraints can be a useful tool for many input design problems and particularly that they can be used to design optimal input signals that are robust with respect to the unknown true system.

Remark: The term variance is used for the asymptotic  ˆ jω ) = lim NE |G(e jω , θˆN ) − G(e jω , θo )|2 . variance Var G(e N→∞

3.2 Input design The input design problem we consider in this paper is the following:  1 π Φu (ω )dω minimize 2π −π (6) jω ˆ subject to Var G(e ) ≤ γ (ω ), ω ∈ Ω.

Before the input design problem is stated, we present some preliminaries regarding system identification and the accuracy of the identified model.

3.1 System identification and model uncertainty

where Ω is some subset of the interval (−π , π ]. The integral of the input spectrum is the total input energy inserted into the system. The objective is to find the input spectrum that minimizes the total input power and still guarantees a certain accuracy of the resulting transfer function estimate.

We will work within the prediction error framework, see e.g (Ljung, 1999). A linear time-invariant (LTI) system is parameterized by the vector θ and the model can be described by the relation y(t) = G(q, θ )u(t) + H(q, θ )e(t)

(2)

The variance constraint in (6), which can be written as γ − Λ∗G PΛG ≥ 0, cf. (5), will be reformulated as a matrix inequality that is linear in P−1 . First view γ −Λ∗G PΛG as the Schur complement of P−1 in 

where y(t) is the output, u(t) is the input, e(t) is white noise with variance λ and q is the forward shift operator, i.e qy(t) = y(t + 1). Further, we assume that the true system can be described within the model structure, represented by a true parameter vector θo and a white noise signal eo (t) with variance λo . We denote by θˆN , the parameter estimate based on N samples of the input/output.

P−1 λG Λ∗G γ

and then consider the Schur complement of γ , which is P−1 − 1 ΛG Λ∗ . Then we see that Var Gˆ ≤ γ is equivalent with

 π

−π

Ro =

1 2π

 π −π

Fe (θo )Fe (θo )∗ dω ,

(7)

In the following sections, the usefulness of using sum of squares constraints in input design problems will be illustrated with a series of examples. The examples will be different modifications and versions of (6). We emphasize that (6) is a prototype example and that the presented methods are easily adapted to some other problem formulations.

Fu (θo )Fu (θo )∗ Φu (ω )dω + Ro (3)

where

G

which is linear in P−1 .

Important here is that P−1 is an affine function of the input spectrum Φu and in the open loop case, P−1 is given by 1 2πλo

γ

P−1 − 1γ ΛG Λ∗G ≥ 0

Provided that enough data is used for the identification, the statistical properties of the model error can be described entirely by the asymptotic covariance matrix of the parameter estimate, here denoted P  limN→∞ NE[(θˆN − θo )(θˆN − θo )T ].

P−1 =

and

3.3 Parametrization of the input spectrum In order to solve the optimization problem (6) we need a parametrization of the spectrum Φu (ω ). In general a spectrum can be expressed as

(4)

1353

Φu (ω ) =

∞

∑

k=−∞

c|k| Bk (e jω )

The basis functions in (8) are Bk = e− jω k and the objective function of the input design problem (6) is simply  1 π Φu (ω )dω = r0 . (14) 2π −π

(8)

B∗−k

where Bk = are some basis function spanning L2 . The notation c|k| is used to emphasize that c−k = ck . There are some different approaches to get a finite dimensional parametrization of the spectrum, see (Jansson, 2004).

4.1 The variance bound

In the first two examples, we use the one referred to as ’partial correlation parametrization’. In that approach, the optimization is performed over the m + 1 first coefficients {ck }m k=0 and to ensure that there exists an infinite expansion {ck }∞ k=m+1 such that jω ) ≥ 0, the Toeplitz matrix c B (e ∑∞ k=−∞ |k| k ⎞ ⎛ c0 c1 · · · cm ⎜ c1 c0 · · · cm−1 ⎟ (9) T =⎜ .. . . . ⎟ ⎝ ... . .. ⎠ . cM cM−1 · · · c0

In this section it will be shown how the variance bound in (6) can be treated as sum of squares constraints. First note that, in this example, Fu (e jω , θ ) = ΛG (e jω , θ ) = Γn(e jω )  [ e− jω ··· e− jω n ]T and that R0 = 0. The inverse of the parameter covariance matrix is given by P−1 =

In the third example we use a finite dimensional spectrum parametrization m

∑

k=−m

c|k| Bk (e jω )

M(e jω )  Ru −

(10)

of the positive real part of the input spectrum, 12 c0 + jω ∑M k=1 ck Bk (e ). The condition Φu ≥ 0 is equivalent with the LMI

    Q−AT QA −AT QB 0 CuT ≥0 K(Q, CAuu DBuu )  −BTuQA u −BuT QBu + C D +D T u u u u u u u (11) T for some

matrix Q = Q . The state space description Au Bu Cu Du can be constructed such that the parameters {ck } appears linearly in Cu and Du . With these approaches, both the objective function and the variance constraint in (6) can be formulated as affine functions of the parameters {ck }m k=0 and (6) can be stated as a finite-dimensional convex optimization problem. More details will be given in the following examples.

∑

k=−∞

rk e− jω k dω ⎞ · · · rn−1 · · · rn−2 ⎟ . ⎟. .. . .. ⎠ · · · r0 (15)

λ0 Γn(e jω )Γn(e jω )∗ ≥ 0 γ

(16)

where the condition T = Ru ≥ 0 is automatically included. M ∈ Cn×n is a complex-valued Hermitian matrix, affine in the real parameters rk , that should be positive semi-definite for all frequencies ω ∈ (−π , π ]. This is equivalent to the larger, real-valued symmetric matrix,     ℜ(M) −ℑ(M) ∈ R2n×2n , (17) M ℑ(M) ℜ(M) being positive semi-definite for all frequencies. ℜ(M) and ℑ(M) means the real and imaginary parts of M. Next, it will be shown how the methods in Section 2 can be used to transform the infinite dimensional constraint  jω ) ≥ 0, ∀ω ∈ (−π , π ] M(e (18) In (18) the unit circle is parameterized by e jω , but it can also be described by the real numbers a and b as a + jb where a2 + b2 − 1 = 0. With e jω replaced by a + jb and e− jω with a − jb, the condition (18) can be written as  + jb) ≥ 0, ∀a, b ∈ R : a2 + b2 − 1 = 0 (19) M(a  now are polynomials in a where the elements of M and b. A sufficient condition for (19) to hold is the existence of a polynomial τ¯ (y, a, b) such that

In this first example we consider input design for identification of an FIR-system. The model is given by

 + τ¯ (y, a, b)(a2 + b2 − 1) = sos(y, a, b) yT My

n

R2n .

(20)

Let the auxiliary polynomial τ¯ (y, a, b) where y ∈ be linearly parameterized by the parameters τk and add them to the free optimization variables. The structure and degree of τ¯ (y, a, b) is quite crucial. If the structure is too limited, the design will be conservative or there may not exist any feasible solutions at all. On the other hand, increasing the degree and flexibility of τ¯ (y, a, b) will increase the dimension of the LMI and also the number of optimization variables will increase.

(12)

k=1

and H(q, θ ) ≡ 1. In this case the spectrum parametrization can be chosen particularly simple. Let ck be the auto-correlation coefficients of the input signal u(t) and in this specific case we denote the parameters rk ck = rk  E{u(t)u(t − k)}.

⎛

∞

into a sum of squares constraint.

4. EXAMPLE 1: FREQUENCY-WISE VARIANCE BOUNDS FOR FIR-SYSTEM

∑ θk q−k

−π

Γn(e jω )Γn(e jω )∗

The variance constraint can be written as, (cf. (5) and (7)),

where make sure that Φu (e jω ) ≥ 0 for all ω .

we must Au Bu Let Cu Du be a controllable state space description

G(q, θ ) =

 π

r0 r 1 ⎜ r1 r0 1 = Ru , where Ru  ⎜ .. ⎝ ... λ0 . rn−1 rn−2

must be constrained to be positive semi-definite. This applies to basis functions Bk (e jω ) = L(e jω )e− jω k where L(e jω ) > 0.

Φu (e jω ) =

1 2πλ0

(13)

1354

The formulation (20) can easily be extended to frequency dependent variance bounds γ (e jω ) by considering bounds on the form |γn (e jω )|2 γ (e jω ) = , |γd (e jω )|2

Bode Diagram 20 15 10

Magnitude (dB)

5

(21)

where γn and γd are polynomials in e jω , and redefining jω

jω

jω

2

jω

2

-20

M(e )  |γn (e )| Ru −|γd (e )| λ0 Γn(e )Γn(e ) . (22)

-25 -30

Define f (a, b)  a2 + b2 − 1 and the input design problem can now be formulated as a sum of squares problem rk ,τk

 + τ¯ (y, a, b) f (a, b) = sos(y, a, b). subject to yT My (23) In this example, let n = 3, λ0 = 1 and let the vari  jω +0.8e2 jω 2  ance bound be γ (e jω ) =  1+0.5e  . Now to j ω 2 j ω 1+0.7e +0.2e the choice of structure of the polynomial τ¯ (y, a, b).  is a quadratic form in y and First notice that yT My contains powers of a and b up to degree 6. To match that in τ¯ (y, a, b) f (a, b) we could choose to include all possible monomials on the form yi y j ak bl , where the product kl ≤ 4, resulting in 540 additional optimization variables. Here we show two less flexible structures, one that gives a conservative design and one that gives the true optimal design. The first one is given by

k=1 l=1

6

yk yl + b4 ∑

(24)

6

(2)

∑ τkl

k=1 l=1

ˆ j ω ) ≤ γi , Var G(e

yk yl

∀ω ∈ (−π , π ].

i = 1, 2, ... (28)

(29)

Below, we will show that these conditions can be handled with sum of squares but also with a generalization of the KYP-lemma (Iwasaki and Hara, 2005).

5.1 Sum of squares The constraints (29) are restricted to segments of the unit circle. Again, consider the unit circle parametrization e jω = a+ jb where f (a, b) = a2 +b2 −1 = 0. Now a segment of the unit circle, ω ≤ ω ≤ ω , can also be described by cos ω ≤ cos ω ≤ cos ω , where cos ω = a and this can be written as one quadratic constraint

For this example, there is an alternative method for solving the optimization problem rk

ω i ≤ ω ≤ ω i,

Mi (e jω ) = γi Ru − λ0 Γn(e jω )Γn(e jω )∗ ≥ 0, ω i ≤ ω ≤ ω i , i = 1, 2, ...

4.3 Reference method via KYP-lemma

subject to M(e jω ) ≥ 0,

3

where the specifications are made in the frequency range ω ∈ (0, π ) and due to symmetry they will also be fulfilled for −ω i ≤ ω ≤ −ω i . The conditions (28) are equivalent to, (cf. (16)),

(25) which has 72 parameters. The results are plotted in Figure 1. The design using (25) is actually the solution to the original problem, i.e. using the constraint (18). Using (24) gives a conservative design where the objective function, (i.e. the input power), is 12% higher than for the optimal design.

minimize r0

2.5

Consider again the FIR-system (12). The only difference in this example from the previous one is the shape of the variance specifications. Here we consider piecewise constant bounds

with only 6 additional parameters and the second one is given by (1)

2

5. EXAMPLE 2: PIECEWISE CONSTANT VARIANCE BOUNDS FOR FIR-SYSTEM

k=1

6

1.5

Frequency (rad/sec)

for some matrix P = PT . Applying this method to our example gives exactly the same solution as using the sum of squares constraint in the previous section, where the auxiliary polynomial τ¯ is given by (25). This comparison is made to certify that the sum of squares method gives a correct result in this case.

6

∑ τkl

1

lemma, see e.g (Yakubovich, 1962), the condition (16) is equivalent to the LMI   P − AT PA −CT C −AT PB ≥0 (27) −BT PA Ru − BT PB

minimize r0

6

0.5

Fig. 1. The solution to the input design problem (23). Solid line – the variance bound γ (e jω ). Dashed line – variance from design using (25). Coincides with the results from Section 4.3. Dotted line – variance from design using (24).

4.2 The input design problem – an example

τ¯ (y, a, b) = a4 ∑

-10 -15

jω ∗

τ¯ (y, a, b) = (a4 + b4 ) ∑ τk y2k

0 -5

(26)

g(a) = ηw2 − (a − ηc )2 ≥ 0

In (Jansson, 2004), similar problems are solved using the KYP-lemma in the following fashion. Express λ0 ∗ ∗ γ ΓnΓn as G G where G has a state space description G(e jω ) = C(e jω I − A)−1 B. Then, by the KYP-

(30)

where ηc = (cos ω + cos ω )/2 and ηw = (cos ω − cos ω )/2. A sum of squares relaxation of (28) can be expressed as

1355

i (a + jb)y + τ¯i (y, a, b) f (a, b) − μ¯ i (y, a, b)g(a) yT M = sos(y, a, b) μ¯ i (y, a, b) = sos(y, a, b), i = 1, 2, ... (31)

Bode Diagram 1

0.9

Magnitude (abs)

0.8

i is defined analogously to (17) and τ¯i and μ¯ i where M are polynomials in y, a and b.

0.7

0.6

0.5

0.4

0.3

5.2 The generalized KYP-lemma

0.2

We will show that the generalized KYP lemma (Iwasaki and Hara, 2005) can be used to obtain nonconservative convex equivalents to the conditions in (28). Express λ0 ΓnΓn∗ as G∗ G where G has a state space description G(e jω ) = C(e jω I − A)−1 B. Then the condition γ Ru − G∗ G ≥ 0 for all ω such that ω ≤ ω ≤ ω is equivalent with the LMI        T −P e j ϑc Q −C C 0 AB ∗ AB ≥0 − − j ϑ c I 0 I 0 0 γ Ru e Q P−(2 cos ϑw )Q (32) for some matrices P = PT and Q = QT ≥ 0 where ϑc = (ω + ω )/2 and ϑw = (ω − ω )/2.

0.1

0.5

1

For the sum of squares approach, the auxiliary polynomials τ¯i and μ¯ i all have the same structure given by k=1

m 1 rk Πk (θo ) ∑ f (θo ) k=−m

(35)

where f and the matrices Πk are polynomial functions of the true parameter vector and where f (θ ) > 0 for all stable θ . Consider a variance bound γ (e jω ) on the form (21) and multiply (7) with |γn |2 |A|4 f . The variance constraint can be formulated as

2

k=1

3

b

Let n = 5, λ0 = 1 and let the variance bound be  γ = 0.5, 0 ≤ ω ≤ 1 γ = γ1 = 0.2, 1 ≤ ω ≤ 2. (33)

(1) (2) τ¯ (y, a, b) = a2 ∑ τk y2k + b2 ∑ τk y2k

2.5

for the AR-process 1/A2 , the inverse of the parameter covariance matrix, cf. (3), can be written

5.3 Numerical example

10

2

Fig. 2. Results from Example 2. Solid line – the variance bound γ (e jω ). Dotted line – variance from design using the sum of squares relaxation (31). Dashed line – variance from optimal design.

1 −BΓna and, solving the Yule-Walker equations AΓn A2

P(θo )−1 =

10

1.5

Frequency (rad/sec)

M(θo ) = |γn |2 |A(θo )|4

m

∑

rk Πk (θo )+ k=−m



−B(θ )Γ −B(θo )Γna ∗ − |γd |2 f (θo )λo A(θoo)Γnna A(θo )Γnb b

(34)

resulting in 80 additional optimization parameters. The resulting variance from this design together with the variance bound is shown in Figure 2. The optimal design computed using the KYP-approach is also shown. We see that in this example the sum of squares relaxation gives a conservative design. The KYP-approach is also significantly faster to compute, so for this type of problem, this approach is to be preferred. However, notice that the objective function obtained with sum of squares is only 2% higher than for the optimal design. Thus, the example serves as an illustration of the potential of sum of squares solutions.

≥ 0, (36)

where the arguments e jω are omitted for brevity. The elements of M(e jω , θo ) are polynomial functions of  j ω , θo ) the parameters θo . Define the real matrix M(e analogously to (17). It is quite clear that the optimal input spectrum will depend on the true unknown system parameter θo , which is an inherent property of most input design problems. If we have some a priori knowledge of the true system, this can be incorporated in the design. Suppose that we know, (or that it is likely), that θo ∈ Θ where Θ is an ellipsoidal region, centered in θ  , given by   (37) Θ = θ : (θ − θ  )T Ω(θ − θ  ) ≤ 1 .

6. EXAMPLE 3: ROBUST INPUT DESIGN FOR OUTPUT-ERROR SYSTEM

Note that (37) is quadratic in θ . Now consider the problem of finding the input with least energy, that guarantees that the variance bound γ is fulfilled for all systems in the set Θ. This can be formulated as

The optimal input spectrum that solves the problem stated in (6) will depend on the unknown true system. We will now show that the sos approach can be used to solve a robustified version of (6) which guarantees that the variance constraint holds for all systems corresponding to an ellipsoidal uncertainty region in the parameter space.

minimize c0 ck ,Q

 jω , θ ) ≥ 0, ∀ω , subject to M(e

K(Q, CAuu DBuu ) ≥ 0.

For simplicity we will consider an output-error system θ) given by G(q, θ ) = B(q, A(q,θ ) and H(q) = 1. In this case we will use the finite dimensional parametrization of the spectrum (10) with Bk (e jω ) = e− jω k , where the input power is given by c0 . Here Fu = ΛG =

∀θ ∈ Θ

(38)

This problem can be solved by a sum of squares relaxation in the following way: Treat the frequency dependance in the same way as in the previous examples. Define h(θ ) = 1 − (θ − θ  )T Ω(θ − θ  ) and consider the optimization problem

1356

techniques can be useful for input design purposes, but much remains to be done. Numerical complexity is certainly a big issue for these methods and that has to be investigated further.

Bode Diagram 20 15 10 5

Magnitude (dB)

0

REFERENCES Bombois, X., G. Scorletti, M. Gevers, R. Hildebrand and P. Van den Hof (2004a). Cheapest openloop identification for control. In: Proceedings 43rd IEEE Conference on Decision and Control. Bahamas. Bombois, X., G. Scorletti, M. Gevers, R. Hildebrand and P. Van den Hof (2004b). Least costly identification experiment for control. In: MTNS’04. Leuven, Belgium. Gevers, M. and L. Ljung (1986). Optimal experiment designs with respect to the intended model application. Automatica 22, 543–554. Goodwin, G.C., J.S. Welsh, A. Feuer and M. Depich (2006). Utilizing prior knowledge in robust optimal experiment design. In: 14th IFAC Symposium on System Identification, SYSID 2006. Hildebrand, R. and M. Gevers (2003). Identification for control: Optimal input design with respect to a worst case ν -gap cost function. SIAM Journal on Control and Optimization 41(5), 1586–1608. Hjalmarsson, H. (2004). From experiment design to closed loop control. Automatica 41, 393–438. Hjalmarsson, H., M. Gevers and F. De Bruyne (1996). For model-based control design, closed-loop identification gives better performance. Automatica 32(12), 1659–1673. Iwasaki, T. and S. Hara (2005). Generalized KYP lemma: Unified frequency domain inequalities with design applications. IEEE Transactions on Automatic Control AC-50(1), 41–59. Jansson, H. (2004). Experiment design with applications in identification for control. Doctoral thesis. KTH. Stockholm, Sweden. Jansson, H. and H. Hjalmarsson (2005). Input design via LMIs admitting frequency-wise model specifications in confidence regions. IEEE Transactions for Automatic Control 50(10), 1534–1549. Ljung, L. (1999). System Identification: Theory for the user. second ed.. Prentice Hall. Löfberg, J. (2004). YALMIP : A toolbox for modeling and optimization in MATLAB. In: Proceedings of the CACSD Conference. Taipei, Taiwan. Available from http://control.ee.ethz.ch/ ˜joloef/yalmip.php. Mårtensson, J., H. Jansson and H. Hjalmarsson (2005). Input design for identification of zeros. In: Proceedings of the 16th IFAC World Congress on Automatic Control. Parrilo, P.A. (2003). Semidefinite programming relaxations for semialgebraic problems. Mathematical Programming 96(2), 293–320. Parrilo, P.A. and B. Sturmfels (2003). Minimizing polynomial functions. Algorithmic and quantitative real algebraic geometry, DIMACS Series in Discrete Mathematics and Theoretical Computer Science 60, 83–99. Yakubovich, V. A. (1962). Solution of certain matrix inequalities occurring in the theory of automatic control. Docl. Acad. Nauk. SSSR pp. 1304–1307. Yuan, Z.D. and L. Ljung (1985). Unprejudiced optimal open loop input design for identification of transfer functions. Automatica 21(6), 697–708. Zhu, Y and P.P.J van den Bosch (2000). Optimal closed-loop identification test design for internal model control. Automatica 36, 1237–1241.

−5 −10 −15 −20 −25 −30 −1

0

10

10 Frequency (rad/sec)

Fig. 3. Results from Example 3. Dotted line – variance bound γ . Solid lines – variance using the robust optimal design for 11 different systems in the uncertainty region minimize c0 ck ,Q,τk ,σk

 + jb, θ )y − τ¯ (y, a, b, θ ) f (a, b)+ subject to yT M(a − σ¯ (y, a, b, θ )h(θ ) = sos(y, a, b, θ ) σ¯ (y, a, b, θ ) = sos(y, a, b, θ )

K(Q, CAuu DBuu ) ≥ 0, (39) where τk and σk are coefficients of the auxiliary polynomials τ¯ and σ¯ . This solution satisfies the variance constraint for all systems in Θ. 6.1 Numerical example In this numerical example we let the model be a simple −1 first order system G(q, θ ) = 1+qθ q−1 , H(q, θ ) = 1. Let

λo = 1 ,γ =

1 |1−0.65e jω |2

and m = 7. Let the uncertainty

region be defined by h(θ ) = 1 − 100(θ + 0.5)2 ≥ 0. In this case, M(e jω , θo ) is a real scalar and the optimization problem we solve is minimize c0 ck ,Q,τk ,σk

subject to M(a + jb, θ ) − τ¯ (a, b, θ ) f (a, b)+ − σ¯ (a, b, θ )h(θ ) = sos(a, b, θ ) σ¯ (a, b, θ ) = sos(a, b, θ )

K(Q, CAuu DBuu ) ≥ 0,

(40)

where the polynomials τ¯ and σ¯ are on the form

σ¯ (a, b, θ ) = τ¯ (θ ) =

7

∑ (σk + (a2 + b2 )σk+7 )θ 2(k−1)

k=1 3

∑ τk θ

(41) 2(k−1)

.

k=1

The result is shown in Figure 3, where the variance using the robust optimal design is plotted for 11 different systems within the uncertainty region, (θ = −0.6 + 0.2k, k = 0, . . . , 10). The variance constraint is satisfied for all 11 systems. 7. CONCLUSIONS We have shown that sum of squares constraints can be used to formulate different kinds of input design problems. This is a first step in evaluating if using these

1357