C IRCUITS S YSTEMS S IGNAL P ROCESSING VOL . 26, N O . 3, 2007, P P. 293–310

c Birkh¨auser Boston (2007)
DOI: 10.1007/s00034-005-0729-z

F REQUENCY D OMAIN A NALYSIS OF U NCERTAIN T IME -VARYING D ISCRETE -T IME S YSTEMS * Przemyslaw Orlowski1

Abstract. The paper develops tools and methods for the analysis of uncertain linear timevarying, discrete-time systems. Our considerations are concerned with the analysis of a system in the frequency domain, in particular with obtaining approximate Bode amplitude diagrams by means of two methods: one based on singular value decomposition and the other on the discrete Fourier transform (SVD-DFT). Starting with the literature review of frequency domain methods for time-varying systems, we introduce the state space model with its equivalent operator form description and the most essential definitions for the Bode diagram approximations using two methods called the spectral and the averaged transfer function, respectively. Next we analyze a system with an additive uncertainty, in particular, given as a linear combination of perturbations. Theoretical considerations are summed up and illustrated by numerical examples of applications to robot manipulators. The obtained bounds (computed by the proposed method) are compared with the results obtained by means of the MATLAB wcgain function. Key words: Uncertain systems, time-varying systems, discrete-time systems, uncertainty analysis, nonstationary systems.

1. Introduction The computation of the frequency bounds of an uncertain system is an essential step in frequency domain based robust control design methods. It is involved in all stages of the design procedure including modeling, system analysis, controller design, and the final design evaluation. Unfortunately, the well-developed concepts and analytic methods of time-invariant systems are not applicable to linear time-variant (LTV) systems. An important issue in robust control is the identification of a set of frequency functions representing the uncertain system, i.e., the set of values of the complex ∗ Received July 29, 2005; revised May 10, 2006; published online June 9, 2007. 1 Szczecin University of Technology, Institute of Control Engineering, Sikorskiego 37, 70-313

Szczecin, Poland. E-mail: [email protected]

294

O RLOWSKI

frequency function at a given frequency. We have developed a method for the identification of value sets based on Bode diagrams. This work was already reported in [18], [17]. It was shown that the proposed Bode diagram sets are efficient ways to describe dynamical systems. They preserve many properties of classical Bode diagrams and are almost the same as those used for LTI systems [18], [17]. Moreover, the basic properties of singular value decomposition (SVD) allow us to eliminate some measurement noise [10], [7]. Some similar frequency domain methods, extensively studied by predecessors, rely on extension of the Laplace transform to varying impulse responses [21], [3], [6], [8]. The result of this analysis is a two-dimensional transfer function (2DTF), which can be further converted into, e.g., an averaged one-dimensional transfer function, as proposed in Section 3.1. Other frequency domain methods are based on the ideas of varying eigenvalues, varying natural frequencies [4], [2], [13], or similar concepts of pseudo-modal parameters [9], [11]. Applications of frequency concepts to uncertain systems have received some attention during past years. Some works are concerned with an analysis of uncertain systems, e.g., Rantzer and Andersson analyze frequency-dependent errors for truncated uncertain LTI models [1]. A probabilistic model in the mediumfrequency range for uncertain systems has been presented by Soize [20]. Techniques for assessing the effect of parametric uncertainties on the higher-order frequency response functions of a nonlinear single degree of freedom system have been studied by Manson and Worden [14]. A comparison of methods for computing the frequency response of uncertain plants was done by Chen and Ballance in [5]. Quite a different approach, concerning mostly robust control system synthesis, can be found in [12]. Other concepts are regions of the complex plane called value sets, which quantify the amount of model uncertainty [19], [16], [15]. This paper introduces a new approach to the analysis of uncertain LTV systems based on the approximated frequency SVD-discrete Fourier transform (DFT) Bode diagrams introduced in [18], [17] and on the averaged method of the 2D TF. This method can be applied to a general uncertain structure including non-linear ones.

2. Model description In order to describe uncertain time-varying discrete-time systems one can use the generalized description employing state equations with time-dependent matrices and additive perturbation matrices. The input-output relationships in real systems are always featured by a nonzero time delay. In such a case the system matrix D(k) ≡ 0 and the term D(k) · u(k) in equation (2) can be omitted.

U NCERTAIN T IME -VARRYING D ISCRETE -T IME S YSTEMS

295

The state equations take the following form: x (k + 1) =(A(k) + A ) · x (k) + (B(k) + B ) · u(k), y (k) =(C(k) + C ) · x (k) + D(k) · u(k), k ∈ N, x(0) = 0,

(1) (2)

where {x, x (k) ∈ Rn , k ∈ {0, . . . , N − 1}} are the nominal and the perturbed state, {u(k) ∈ Rm , k ∈ {0, . . . , N − 1}} is the nominal control, {y, y (k) ∈ R p , k ∈ {0, . . . , N − 1}} are the nominal and the perturbed output, and {A(k), A ∈ Rn×n , B(k), B ∈ Rn×m , C(k), C ∈ R p×n , k ∈ {0, . . . , N − 1}} are system matrices with perturbations. The system (1), (2) is called perturbed if at least one of the perturbation matrices is nonzero. If all perturbation matrices are identical to zero the system is called nominal. The output and state trajectories of the perturbed system are marked by the subscript , which allows us to distinguish between the nominal (without a subscript) and the perturbed () trajectory. The model of the system can be described by means of the following operators:   0 0 ··· 0 0  I 0 ··· 0 0    .. ..   ˆ A(1) I 0 . . L= ,   . . . .  . . I 0 0 (3) A(N − 2) · . . . · A(1) · · · A(N − 2) I 0   I   A(0)  ˆ = N   ..   . A(N − 2) · . . . · A(0)     B(0) 0 0 C(0) 0 0    .. .. ˆ = Bˆ =  0 C (4) ,  0 , . . 0 0 0 0 B(N − 1) 0 0 C(N − 1) ˆ have block diagonal forms. The states xˆ , xˆ  , outputs ˆ and C where the operators B ˆ Lˆ B ˆ is yˆ , yˆ  , and the input uˆ have block column vector forms. The operator C a compact, Hilbert-Schmidt operator from l2 into l2 and actually maps bounded signals u(k) ∈ U = l2 [0, N ] into bounded signals y ∈ Y. The output trajectory of the nominal system can be described in the following form: ˆ ·N ˆ ·L ˆ · x0 + C ˆ ·B ˆ · u. ˆ yˆ = C

(5)

It is assumed in this paper that the initial conditions are equal to zero, and thus the system response with zero, initial conditions is determined by the term yˆ u = ˆ ·L ˆ · Bˆ · u. ˆ C

296

O RLOWSKI

3. Bode characteristics approximation 3.1. Averaged transfer function method One of the first attempts of analysis of LTV systems in the frequency domain was made in [21] and developed by followers, e.g., [3], [6], [8]. The time-varying transfer function has been defined by extending the Laplace transform to the varying impulse response. The most general realization of continuous-time systems is the following generalized Weyl symbol (GWS) introduced by Kozek: 





H (α) (t, f ) = h t + 12 − α τ, t − 12 + α τ e− j2π f τ dτ, (6) τ

where α ∈ R is an arbitrary real number, usually bounded such that |α| ≤ 0.5 and h(t1 , t0 ) is the response of the system to the Dirac impulse (or Kronecker delta for discrete-time systems) δ(t − t0 ) at time t1 (after t1 −t0 seconds). The time-varying transfer function [21] is obtained for α = 0.5, frequency-dependent modulation function [3] for α = −0.5, and the Weyl symbol for α = 0. In general, the GWS does not satisfy a few of the properties that hold for LTI systems. First of all, the cascade connection of two systems cannot correspond to multiplication of the components in the frequency domain, which holds for LTI systems (H12 = H1 H2 ⇔ H12 ( f ) = H1 ( f )H2 ( f )). Moreover, in general the input-output relationship does not take the form of a multiplication in the frequency domain (for LTI systems only y(t) = (Hx)(t) ⇔ Y ( f ) = H ( f )X ( f )). For numerical computations equation (6) is discretized and the DFT is used, H (α) (tl , f k ) =

N





 h l + 12 − α n T p , l − 12 + α n T p e− j2π(k−1)(n−1)/N, n=1

(7) where tl = lT p , f k = (k−1) N T p , l = 1, 2, . . . , N , and k = 1, 2, . . . , N /2. It is difficult to analyze properties of the LTV system by directly using the GWS, because the data represent complex surfaces in a 3D space. The simplifications rely on the transformation of these surfaces into 2D diagrams and application of tools similar to those for LTI systems. Definition 1. The averaged transfer function (ATF) for a linear time-varying discrete-time system on the finite time horizon is given as a mean of the discrete 2D transfer function H with respect to time, G (α) A ( jωk ) =

N 1  H (α) (ti , f k ). N i=1

(8) (−0.5)

If α is unspecified it is assumed that α = −0.5 and G A ( jωk ) = G A ωk = 2π f k .

( jωk ),

U NCERTAIN T IME -VARRYING D ISCRETE -T IME S YSTEMS

297

The definition follows from the property satisfied for LTI systems, H (α) (tl , f k ) = H (α) (tm , f k ) =

N 1  H (α) (ti , f k ), N i=1

(9)

where l and m are arbitrary integers in a given range l, m ∈ {1, 2, . . . , N }. The definition of ATF for an LTV system relies upon the assumption that every sample has the same importance and thus can be averaged. An important property of the ATF is linearity (additivity and homogeneity). It arises from linearity of both the DFT and the mean. The ATF can be divided (α) into the magnitude and the phase plot, where the magnitude |G A ( jωk )| can be interpreted as the selective amplification of the first harmonic in the output spectra for a given sinusoidal input. All the other components which exist only for timevarying systems are neglected. 3.2. Spectral transfer function method Some elements of the frequency analysis based on SVD-DFT decomposition of the system operator have been introduced in [18]. The general result of the analysis is the spectral transfer function (STF), which is recalled in the separate magnitude and phase forms. The magnitude-frequency response |G(ωk )| is uniquely defined by

N

1  σ 2 · |DFTk [u j ]|2 |G S (ωk )| =  · N j=1 j (10)  2

N  N 

1     = · u ni · e− j·2·π·(k−1)·(n−1)/N  , σi ·  N i=1  n=1 and the phase-frequency response ϕ(ωk ) = arg(G(ωk )) can be written as   N  DFTk [u j ] ϕ S (ωk ) = arg σj · DFTk [v j ] j=1 N  N  − j·2·π·(k−1)·(n−1)/N  n=1 u ni · e σi ·  N , = arg − j·2·π·(k−1)·(n−1)/N i=1 n=1 vni · e

(11)

where U · S · VT is the product of SVD of the system input-output operator ˆL ˆ B, ˆ such that S = diag(σi ) is a diagonal matrix with singular values and U, V C are orthonormal matrices composed of column vectors ui and vi , respectively. Singular values σi in equations (10), (11) play the role of weight functions. The derived relationships hold true for both the time-invariant and the time-variant systems. The classical Bode diagrams calculated for LTI systems on the infinite time horizon (ITH), using the substitution z = exp( jωT p ) in the discrete transfer

298

O RLOWSKI

function, are almost identical for either the ATF or STF method. Affinity of the diagrams holds only for a large enough finite time horizon (FTH) used in ATF/STF methods. In contrast to ATF, STF is only homogeneous but not additive. This follows from equation (10) and the quadratic mean, which is nonlinear. The STF can also be divided into the magnitude and phase plot, but in this case the magnitude |G S ( jωk )| can be interpreted as a cumulative amplification of all harmonics in the output spectra for a given sinusoidal input. Despite the nonlinearity, this method gives better results, e.g., for stability analysis of the feedback loop.

4. Uncertainty analysis In order to analyze uncertain LTV systems, one has to introduce a few new opˆ A,  ˆ B,  ˆ C that have block diagonal forms similar to equation (4). Iˆ erators,  denotes the identity operator (matrix). Using this operator notation, one can prove the following theorem for uncertain discrete LTV systems. Theorem 1. For any uncertain discrete-time LTV system described by the nomiˆL ˆB ˆ and the perturbation additive operators given in nal input-output operator C the following form: ˆ = A ˆ + ˆ = C ˆ + ˆ = B ˆ + ˆ A, ˆ B, ˆ C, A B C (12) ˆ Lˆ B) ˆ  can be written in terms the input-output operator of the perturbed system (C of the nominal operator and the additive perturbation operator as follows: ˆL ˆL ˆB ˆ + ˆ B) ˆ =C ˆ CLB , (C (13) where the perturbation operator is given by ˆ  (L ˆ  Lˆ  ˆ ˆ ˆ CLB = C ˆB+ ˆ A + (L ˆ A )2 + · · · )Lˆ Bˆ  + C ˆ C Lˆ Bˆ 

(14)

and the following condition holds: ˆ · ˆ A  < 1. L

(15)

Proof. The nominal state trajectory in the case of zero initial conditions can be written as follows: ˆ ·B ˆ · u. ˆ xˆ p = L (16) It has been proven [16] that the perturbed state trajectory can be written as ˆ · ˆ · ˆ A · xˆ  + L ˆ B · u. ˆ xˆ  = xˆ p + L

(17)

ˆ · ˆ A ) has to be invertible. The necessary condition for The expression (Iˆ − L ˆ · ˆ A  < 1. Thus, equation (17) can be rearranged into the its invertibility is L following form: ˆ · ˆ ·B ˆ · uˆ + L ˆ · ˆ A )−1 · (L ˆ B · u), ˆ xˆ  = (Iˆ − L (18)

U NCERTAIN T IME -VARRYING D ISCRETE -T IME S YSTEMS

299

and the perturbed output trajectory is given by ˆ + ˆ · ˆ · (B ˆ + ˆ C ) · (Iˆ − L ˆ A )−1 · L ˆ B ) · u. ˆ yˆ  = (C

(19)

ˆ · ˆ A  < 1 it is possible to rewrite the inverted term Under the assumption L ˆ ˆ ˆ · ˆ · ˆ A )−1 = Iˆ + L ˆ A + (L ˆ A )2 + · · · using a Taylor series (I − L ·  ˆ + ˆ · ˆ · ˆ · (B ˆ + ˆ C ) · (Iˆ + L ˆ A + (L ˆ A )2 + · · · ) · L ˆ B ) · u, ˆ yˆ  = (C

(20)

and the output perturbed trajectory can also be written in terms of perturbed operators ˆ L ˆ Lˆ B) ˆ  · u. ˆ B ˆ  · uˆ = (C ˆ yˆ  = C

(21)

Thus, the perturbed operator is given by ˆL ˆL ˆ L ˆL ˆB ˆ + C( ˆ ˆ ˆ Bˆ + C ˆ ˆ B) ˆ  =C ˆ A + (L ˆ A )2 + · · · )L ˆB+ ˆ C Lˆ Bˆ (C ˆ Lˆ  ˆ ˆ ˆ ˆ ˆ A )2 + · · · )L ˆB+ ˆ C (L ˆ A + (L ˆ A )2 + · · · ) ˆ A + (L + C( ˆ ˆ ˆ ˆ ˆ CL ˆ C (L ˆ A + (L ˆ A )2 + · · · )L ˆ B, ˆB+ · Lˆ Bˆ +  (22) and, after simplification, ˆL ˆL ˆ  (L ˆ  Lˆ  ˆ B) ˆ =C ˆB ˆ +C ˆ ˆ ˆ Bˆ  + C ˆ ˆ A + (L ˆ A )2 + · · · )L ˆ C Lˆ B, ˆB+ (C (23) ✷

which is equivalent to (14) and thus finishes the proof.

For numerical computations and small perturbations, the results of Theorem 1 can be simplified in the following way. ˆ ) ˆ n  Corollary 1. For suficiently small perturbation matrices, we have (L 1 for n ≥ 2, and hence all powers of these matrices higher than or equal to 2 can be omitted. By restricting the perturbation matrices to the first-order, the firstorder approximation of the perturbation operator is given by ˆL ˆL ˆ ˆB ˆ +C ˆ ˆ ˆ AL ˆ C Lˆ B. ˆ (1) = C ˆB+  CLB

(24)

Further simplifications are possible under more specific assumptions on the uncertainty matrices. Furthermore, Theorem 1 can be extended to a more general class of uncertain systems called systems with parameter bounded uncertainty of a linear combination. Such a generalization is presented below. Corollary 2. A generalization of the additive uncertainty given in (24) is the uncertainty of a linear combination defined for each perturbation as follows:  = d1 · D1 + · · · + dk · Dk =

k  i=1

di i .

(25)

300

O RLOWSKI

k k k ˆ Ai ,  ˆ Bi ,  ˆ Ci , ˆ A = i=1 ˆ B = i=1 ˆ C = i=1 After substituting  ai  bi  ci  equation (24) can be rewritten in the form ˆ (1) =  CLB

k  i=1

ˆL ˆ ˆB ˆ+ ˆ Ai L ai C

k 

ˆ Lˆ  ˆ Bi + bi C

i=1

k 

ˆ ˆ Ci Lˆ B. ci 

(26)

i=1

ˆ CLB , e.g., For a given number k of terms, one can easily write the relation for  for k = 2 the perturbation operator is given by ˆ + c1  ˆ ˆ CLB =(C ˆ C1 + c2  ˆ C2 )(a1 L ˆ A1 + a2 Lˆ  ˆ A2  ˆ Bˆ + b1  ˆ ˆ ˆ B1 + b2  ˆ B2 ) ˆ A1 + a2 L ˆ A2 )2 + · · · )L( + (a1 L ˆ + c1  ˆ 1 ˆ B. ˆ ˆ C1 + c2  ˆ C2 )L(b ˆ B1 + b2  ˆ B2 ) + (c1  ˆ C1 + c2  ˆ C2 )L + (C (27) For computational reasons the simplified relation (26) is used rather than the relation (27). 4.1. Frequency domain transformation ˆL ˆ Bˆ can be transformed into the frequency domain using Every system operator C the methods described in Section 3. Such a transformation can also be carried out also for the perturbed operator. For LTI systems, the following relation holds: ˆL ˆL ˆ B) ˆ + F( ˆ B) ˆ  ) = F(C ˆ CLB ), F((C

(28)

ˆL ˆL ˆ B)(ω) ˆ ˆ B) ˆ  (ω) = (C ˆ CLB )(ω), + ( (C

(29)

or equivalently, where F can be either the F A or F S transformation defined by equations (8) or ˆL ˆL ˆ B) ˆ = (C ˆ B)(ω), ˆ (10), (11), respectively, and F(C etc. Equation (28) is satisfied for both F A and F S transforms of LTI systems and for the F A transform of LTV systems. In general, it is not satisfied for the F S transform of LTV systems because of the square mean in equation (10) and the properties of the SVD. However, the following relation holds for the F S transform of an LTV system: ˆL ˆL ˆ B) ˆ  )| ≤ |FS (C ˆ B)| ˆ + |FS ( ˆ CLB )|, |FS ((C

(30)

and the following condition always holds: lim

ˆ CLB →0 

ˆL ˆL ˆ B) ˆ  ) − F(C ˆ B) ˆ − F( ˆ CLB )] = 0. [F((C

(31)

Thus, equation (28) is also approximately satisfied for the F S transform of LTV ˆ CLB is sufficiently small. What do we mean by “sufficiently”? systems when  ˆ Lˆ B) ˆ It was tested by numerical computations that for the case F(C = 4· ˆ F(CLB ) under the assumption that the perturbation is time invariant, the

U NCERTAIN T IME -VARRYING D ISCRETE -T IME S YSTEMS

301

error in (28) is approximately equal to 6%. This error decreases when the factor ˆL ˆL ˆ B) ˆ decreases. For both C ˆ B, ˆ  ˆ CLB )/F(C ˆ CLB time varying, δF = F( the error does not exceed the factor δF . 4.2. Uncertainty estimates The deviation of magnitude and phase diagrams can be estimated using a method that is based on simple geometrical relationships for complex numbers. It can be summarized as the following theorem. Theorem 2. Lower and upper bounds of the magnitude and phase plots can be estimated in the following way: ˆL ˆL ˆ Lˆ B) ˆ B) ˆ −1 (ω)| ≤ |(C ˆ B) ˆ  (ω)| ≤ |(C ˆ +1 (ω)| |(C

(32)

ˆ Lˆ B) ˆL ˆ Lˆ B) ˆ −i (ω)) ≤ arg((C ˆ B) ˆ  (ω)) ≤ arg((C ˆ +i (ω)), arg((C

(33)

where ˆL ˆ Lˆ  ˆ ˆ B)(ω)), ˆ ˆ Ai L ˆ Bi )(ω)), v Ai (ω) = arg((C v Bi (ω) = arg((C ˆ B)(ω)) ˆ ˆ Ci L vCi (ω) = arg((   ˆ Lˆ  ˆ ˆ Ai Lˆ B)(ω) a¯ i sgn(cos(v Ai (ω) − vq (ω)))(C k   ¯ ˆ (1) )(ω) = ˆ Lˆ  ˆ Bi )(ω)  (  +bi sgn(cos(v Bi (ω) − vq (ω)))(C CLBq i=1 +c¯ sgn(cos(v (ω) − v (ω)))( ˆ ˆ Ci Lˆ B)(ω) i Ci q

(34)

(35)

ˆL ˆL ˆ B)(ω) ˆ ˆ B) ˆ q (ω) = (C ˆ CLBq )(ω) (C + (

(36)

ˆL ˆ B) ˆ  (ω)) vq (ω) = arg(q · (C

(37)

ˆL ˆ B) ˆ q (ω)  = 0 for all ωk and q = {+1, −1, +i, −i}, (C k = 1, 2, . . . , N /2, and q.

=

2π(k−1) N Tp ,

Equations (32) and (33) can be easily proven geometrically for any fixed ω, for which the DFT can be understood as a complex number. The magnitude of an uncertain system is bounded by the parallel perturbation estimate given by equation (35) with q = {−1, +1}. The estimate is the worst possible combination of ˆL ˆ B) ˆ q (ω) (for −1 equation (26) calculated in the direction parallel to operator (C in the reverse direction). The phase of an uncertain system can be bounded by the orthogonal estimate given by (35) with q = {−i, +i}. The estimate corresponds ˆ Lˆ B) ˆ q (ω). to directions orthogonal to the operator (C ˆ ˆ ˆ If (CLB)q (ω) = 0, then Theorem 2 cannot be used for estimation of the bounds. From a practical point of view, such a singular case can be converted to a ˆL ˆ B) ˆ q (ω)  = 0, e.g., by changing the sampling period or the nonsingular case (C simulation horizon.

302

O RLOWSKI

ˆ A,  ˆ B,  ˆ C the perturbation operator  ˆ CLBq can be Corollary 3. For small  replaced by its first-order approximation, thus

ˆ Lˆ B) ˆL ˆ q (ω) ∼ ˆ B)(ω) ˆ ˆ (1) )(ω) (C + ( = (C CLBq

(38)

ˆL ˆ B) ˆ  in equation (37) can be approximated by the operator and the operator (C ˆL ˆB ˆ and vq is approximated by C ˆL ˆ B)(ω)) ˆ vq (ω) ∼ and q = {+1, −1, +i, −i}. = arg(q · (C

(39)

ˆ CLBq Alternatively, the first-order approximation of the perturbation operator  can be given in the following form:    ˆˆ  ˆ B)(ω) ˆ ˆ Ai L (CL ˆ Lˆ  ˆ ˆ Ai Lˆ B)(ω) (C a¯ i sgn Re q ˆ ˆ ˆ  (CLB)q (ω)    ˆˆ  k   ˆ Bi )(ω)  (CL (1)  ˆ Lˆ  ¯ ˆ Bi )(ω) CLBq (ω) = (C +bi sgn Re q ˆ ˆ ˆ (CLB)q (ω) i=1       ˆ B)(ω) ˆ ˆ Ci L  ( ˆ ˆ Ci Lˆ B)(ω) +c¯i sgn Re q ( ˆL ˆ B) ˆ q (ω) (C

      .      (40)

ˆL ˆ B) ˆ q (ω)  = 0. The form is equivalent to equation (35) under the assumption (C

5. Numerical examples The purpose of the examples is twofold. First, they test the proposed algorithm under different conditions. Second, they illustrate the use of the method. The first example compares the results for the ATF and STF methods with bounds of uncertain systems computed by means of Robust Control Toolbox v.3.0 (MATLAB 7.0.1). The second example shows the superiority of the ATF/STF methods applied to an LTV system on an FTH.

U NCERTAIN T IME -VARRYING D ISCRETE -T IME S YSTEMS

303

5.1. LTI system The system is characterized by the following matrices of the discrete-time statespace model (1), (2):      T 0.2574 0.0234 −0.95 0.57 A(k) = , B(k) = , C(k) = 0.0234 0.2132 0.78 −0.82 (41)     0.0515 0.0047 −0.095 A (k) = α · , B (k) = β · , 0.0047 0.0426 0.078  T 0.28 C (k) = γ · . (42) −0.41 The parameters of the system are as follows: sampling period T p = 0.004 s, time horizon N = 100 steps, α, β, γ ∈ [−1, 1]. The eigenvalues of the matrix A(k) are approximately equal to 0.2, 0.27. The system has one zero approximately equal to 0.26. A comparative analysis has been carried out for the proposed method and the wcgain function from MATLAB Robust Control Toolbox v.3.0. Figure 1 shows Bode diagrams for both methods. The solid thick line shows the diagram calculated for the nominal system. The two thin dashed lines are ATF estimates calculated from equation (32), (33) (for the LTI system ATF is equivalent to STF). The thin dotted lines show diagrams calculated for a few real cases, not ˆ ˆ A )2  ∼ estimates. The calculated norm of the omitted square term (L = 0.0027. The maximal worst-case gain readout from the above estimates for the FTH case is equal to 7.86 dB. In contrast, the worst-case gain computed for the ITH case using wcgain implemented in MATLAB gives two values: 8.31 dB, the lower bound of the worst-case gain, corresponding to the dash-dot line on Figure 1, and 8.67 dB, the upper bound of the worst-case gain, both marked by circles. One can notice a small difference of 0.45 dB, or 5% in the linear scale, between the FTH and ITH (wcgain) solutions. It is caused by the FTH (100 steps), finite resolution in the frequency domain, and some inaccuracies of the SVD-DFT based diagrams (ATF/STF) for frequencies close to zero. The perturbation matrices for the worstcase gain are identical for both cases (α = 1, β = 1, γ = 1). 5.2. Robot manipulator In contrast to the previous subsection, the LTI system will be replaced now by an LTV planar robot manipulator with varying inertia links. The manipulator is placed in a horizontal plan, and each link has a sliding mass whose position can be varied. Such a manipulator was discussed and analyzed using the pseudo-modal parameters method in [11]. The model description and all parameters are gathered from [11] and discretized with period T p = 0.004 s. The system is of the fourth

304

O RLOWSKI

Figure 1. Magnitude- and phase-frequency responses with upper and lower bounds using ATF method and computed values from wcgain function obtained for LTI system with a uniform distribution.

order; the simulation horizon is equal to N = 500. The distinction from previous works is that the system is analyzed not only with nominal values of parameters, but also with small perturbations of the system matrices. The method proposed in the previous section is founded on the uncertainty model given by a linear combination. It would be interesting to analyze the error of the estimates if the uncertainty model is not strictly satisfied. Thus, the model of uncertainty assumed for numerical examples is only approximately equivalent to the model from equation (26). The system has been analyzed for two different cases of perturbation matrices. First, when the perturbation matrices have normal distribution N , A (k) = {δi,A j (k) = σ A · ai, j (k) · N (0, 1)} or

(43) A (k) =

{ai,j (k)

= N (ai, j (k), (σ A · ai, j (k)) )} 2

and second, when the perturbation matrices have rectangular distribution , A (k) = {δi,A j (k) = σ A · ai, j (k) · (−1, 1)} or

(44) A (k) =

{ai,j (k)

= ai, j (k) + σ A · ai, j (k) · (−1, 1)},

U NCERTAIN T IME -VARRYING D ISCRETE -T IME S YSTEMS

305

where A = {ai, j }, N (µ, σ 2 ) is the normal distribution with mean µ and variance σ 2 , and (−1, 1) is the uniform distribution in the range −1, 1. Similar relationships hold for other system matrices B, C. For the purpose of the examples, the following notions are used. The first one is the notion of timevarying perturbation (TVP) estimates. The perturbation matrices A , B , C are time varying and calculated for a realization of the system assuming a random distribution of system parameters. For example, the perturbation operator of matrix A is equal to   0 A (0) 0  .. ˆA = A (k) = A (k) − A(k) (45)  ,  0 . 0 0 0 A (N − 1) ˆ B and  ˆ C . Such a form does not and similarly for the perturbation operators  guarantee any specific properties for our estimates but is quite useful for comparison with the TIP estimates described below and also for evaluation of the proposed estimates. The second one is the notion of time-invariant perturbation (TIP) estimates. The perturbation matrices are estimated by averaging system matrices for a realization of the system. The bounds for matrices a¯ i , b¯i , c¯i are taken from the distribution parameters of perturbations, which can be easily estimated in real case situations. For comparison, the perturbation operator of the matrix A is equal to   0 Am 0 N 1   . ˆ Ati p =  .  ,  = σ A(k) (46)  0 A Am . 0  N k=1 0 0 Am ˆ B and  ˆ C . Also, this form does not guarantee any and similarly for the matrices  specific properties for the estimates, but in practice it is a good approximation of real deviations. The amplitude and phase diagrams calculated for the STF, using the method of Section 4 for parameters perturbed with a normal distribution, are depicted in Figure 2. The solid thick line shows a diagram calculated for the nominal system. The two solid thin curves above and below are TVP estimates. The two thin dashed lines are TIP estimates. Both TVP and TIP estimates have been calculated from equations (32), (33). The thin dotted line shows diagrams calculated for a few real cases, not estimates. The parameters assumed in simulations are as follows: σ A = 0.005, σ B = 0.1, σC = 0.02. The calculated norm of the omitted ˆ ˆ A )2  ∼ square term is (L = 0.17. As can be seen in Figure 2, all diagrams calculated for randomly perturbed real systems are close to the nominal one and lie inside the calculated bounds both for TVP and averaged TIP perturbations. Figure 2 shows that the phase temporarily goes outside the estimated bounds. It is caused by significant simplifications for both the model of the uncertainty and the estimation method in the direction associated with equation (39).

306

O RLOWSKI

Figure 2. Magnitude- and phase-frequency responses with upper and lower bounds for LTV system with a uniform distribution obtained using STF method.

Figure 3 also shows approximated Bode diagrams, but in contrast to Figure 2, they have been calculated using the ATF method. The parameters are perturbed with a uniform distribution: σ A = 0.005, σ B = 0.1, σC = 0.1. The calculated ˆ ˆ A )2  ∼ norm of the omitted square term is (L = 0.48. As can be seen, these diagrams are closer to the nominal diagram even in the case of higher values of the perturbation. The estimates for both the STF and ATF methods are outside the real plots depicted with the thick dotted line, even in spite of the assumed model of uncertainty, which is not strictly a linear combination of given perturbations. An exception is the phase plot in Figure 2, for which the relations are not strictly satisfied. Now let us try to compare the results with wcgain. The wcgain in MATLAB can be computed only for LTI systems. Thus, some data in the model must be neglected. It is assumed that the system matrices are invariant (frozen at time k = 0). The analysis can be carried out in two ways.

U NCERTAIN T IME -VARRYING D ISCRETE -T IME S YSTEMS

307

Figure 3. Magnitude- and phase-frequency responses with upper and lower bounds for LTV system with a normal distribution obtained using ATF method.

The first way is to convert variability in the matrices into the uncertainty of the matrices, append the real uncertainty, and then apply the wcgain function. Both the time variability and the uncertainty must be introduced by the ureal function (uncertain real parameters) instead of udyn (uncertain dynamics), which is not currently handled by wcgain. The results are not encouraging; both lower and upper worst-case gains are equal to infinity. Although such results make sense for the ITH case, they are meaningless for the FTH one. The response of the DT-FTH system defined by the Hilbert-Schmidt operator is always bounded, and the gain must be finite. The second possibility is to neglect the variability and apply wcgain similarly as for the LTI systems. By doing so one can compute the lower worst-case gain as 44 dB and the upper equal to infinity. Figure 4 shows the magnitude diagram determined using the ATF and wcgain method. Already the lower worst-case gain gives results much higher than the proposed method. It follows from the fact that

308

O RLOWSKI

Figure 4. Magnitude with upper and lower bounds using ATF method and FTH diagram for the worstcase solution from wcgain obtained for LTV system with a normal distribution.

wcgain uses an ITH, while our method makes computations on an FTH (in this case 100 steps). The two thin dashed lines are ATF estimates calculated from equations (32), (33). The thick dash-dotted line shows the worst case calculated from wcgain. It is significant to notice the difference between the computed value (44 dB) and the real maximal value readout from Figure 4 (−3.4 dB), which is smaller than the ATF estimate (−0.8 dB). This follows from the FTH and the finite resolution in the frequency domain.

6. Conclusion The study has achieved two tasks: 1. Extension of the ATF and STF methods to the analysis of uncertain systems 2. Verification of theoretical results by numerical computations. For the first task the discrete-time state-space model with additive perturbations has been used to represent an LTV uncertain system. The two particular methods for LTV Bode diagrams approximation can be extended to derive estimates for uncertain LTV systems. The frequency bounds for LTV systems can be estimated using the method proposed in Section 4.

U NCERTAIN T IME -VARRYING D ISCRETE -T IME S YSTEMS

309

For the second task some examples have been analyzed using the MATLAB environment. The obtained results show that the calculated estimates are close to the real perturbed diagrams. The main advantage of the proposed method with respect to the wcgain analysis follows from the fact that the effect of uncertainty can be analyzed on a finite, arbitrary time horizon as opposed to an ITH. Moreover, the proposed method computes not only the maximal gain (real scalar) but also the lower and upper bounds for all frequencies (two complex vectors). Although it is not always needed for control synthesis, it is useful for other applications such as the analysis of communication channels where not only the worst case gain is needed. Another important advantage of the proposed method is a new cognitive value introduced into the analysis of uncertain LTV systems. One can not only estimate amplitude and phase diagrams, equivalents of the well-known Bode diagrams, but can also calculate the range of possible realizations of the diagram under the influence of uncertainty existing in the model. The examples from Section 5 show the usefulness of the proposed method and a good convergence of the solution. The solution is convergent not only under simplifications from Section 4 but also in the presence of inaccuracies of the model for simulations and estimations. However, we note that the method also has some imperfections. The first is related to the limited amount of data that can be placed on a 2D plot. This follows from the accepted a priori assumption of simplification and reduction of the system. On the other hand, the analysis of data on 3D plots is much more complicated, and such a simplification is an important advantage of the method in some cases. The second limitation is related to the convergence of the matrix Taylor ˆ · A  << 1) and the assumption (Lˆ · A )n ≈ 0 for n = 2, 3, . . . , series (L which allows us to neglect higher-order terms. From a practical point of view, it ˆ means that the maximal possible perturbation of  A is bounded by operator L and the system can be analyzed only for sufficiently small  A . References [1] L. Andersson and A. Rantzer, Frequency-dependent error bounds for uncertain linear models, IEEE Trans. Automat. Control, 44:11, 2094–2098, 1999. [2] M. Basseville, A. Benveniste, G. Moustakides, and A. Roug´ee, Detection and diagnosis of changes in the eigenstructure of nonstationary multivariable systems, Automatica, 23, 479–489, 1987. [3] P. A. Bello, Characterisation of randomly time-variant linear channels, IEEE Trans. Comm. Syst., 11, 360–393, 1963. [4] N. N. Bogoliubov and Y. A. Mitropolsky, Asymptotic Methods in the Theory of Non-linear Oscillations, Hindustan Publishing Corp., Delhi, 1961. [5] W. H. Chen and D. J. Ballance, A comparison of methods for computing frequency response of uncertain plants, in IFAC World Congress, 30, pages 7–12, 1999. [6] H. D’Angelo, Linear Time-Varying Systems, Allyn & Bacon, Boston, 1970. [7] S. Hashemi and J. K. Hammond, The interpretation of singular values in the inversion of minimum and non-minimum phase systems, Mechanical Systems and Signal Processing, vol. 10, 225–240. 1996.

310

O RLOWSKI

[8] W. Kozek, On the generalized transfer function calculus for underspread LTV channels, IEEE Trans. Signal Proc., 45, 219–223, 1997. [9] K. Liu, Modal parameter estimation using the state space method, Journal of Sound and Vibration, 197, 387–402, 1996. [10] K. Liu, Identification of linear time-varying systems, Journal of Sound and Vibration, 204, 487–505, 1997. [11] K. Liu, Extension of modal analysis to linear time-varying systems, Journal of Sound and Vibration, 226, 149–167, 1999. [12] H. Logemann and S. Townley, Low gain control of uncertain regular linear systems, SIAM J. Control and Optimization, 35, 78–116, 1997. [13] N. M. M. Maia and J. M. M. Silva, Theoretical and Experimental Modal Analysis, John Wiley & Sons, New York, 1997. [14] G. Manson and K. Worden, Frequency response functions for uncertain nonlinear systems, Materials Science Forum, vol. 440–441, 37–44, 2003. [15] M. Milanese, J. Norton, H. Piet-Lahanier, and E. Walter (Eds.), Bounding Approaches to System Identification, Plenum, New York, 1996. [16] P. Orlowski, Deviations estimates for uncertain time-varying discrete-time systems, Third IFAC Symposium on Robust Control Design, Prague, Proceedings art. no. 132, 2000. [17] P. Orlowski, An introduction to SVD-DFT frequency analysis for time-varying systems, Proc. of 9th IEEE Internat. Conf. MMAR, Midzyzdroje, Poland, pp. 455–460, 2003. [18] P. Orlowski, Selected problems of frequency analysis for time-varying discrete-time systems using singular value decomposition and discrete Fourier transform, Journal of Sound and Vibration, vol. 278, issues 4–5, pp. 903–921, 2004. [19] H. Rotstein, N. Galperin, and P. Gutman, Set membership approach for reducing value sets in the frequency domain, IEEE Trans. Automat. Control, 43, 1346–1350, 1998. [20] C. Soize, Uncertain dynamical systems in the medium-frequency range, Journal of Engineering Mechanics-ASCE, 129(9), 1017–1027, 2003. [21] L. A. Zadeh, Frequency analysis of variable networks, Proceedings of the Institute of Radio Engineers, 38, 291–299, 1950.