# Lyapunov Exponents for Bilinear Systems

International Journal of Bifurcation and Chaos, Vol. 13, No. 3 (2003) 713–722 c World Scienti?c Publishing Company

LYAPUNOV EXPONENTS FOR BILINEAR SYSTEMS. APPLICATION TO THE BUCK CONVERTER

CARLES BATLLE? , IMMA MASSANA? and ALICIA MIRALLES? Departament de Matem` atica Aplicada IV, EUPVG and EPSC, Technical University of Catalonia, Av. V. Balaguer s/n, 08800 Vilanova i la Geltr? u, Spain ? carles@mat.upc.es ? imma@mat.upc.es ? almirall@mat.upc.es Received July 12, 2001; Revised February 13, 2002

Using the QR decomposition of the tangent map, we explicitly write the di?erential equations satis?ed by the Lyapunov exponents of a two-dimensional bilinear system, and specialize them to the buck converter in continuous current mode. These equations are solved numerically for an arbitrary value of the bifurcation parameter, and analytically for periodic orbits, recovering for the later results already known from the computation of Floquet exponents. Keywords : Lyapunov exponents; bilinear systems; buck converter.

1. Introduction

Lyapunov exponents measure the average rate of stretching of in?nitesimal elements in state-space, with positive values corresponding to separation of trajectories in a given direction. In?nite-time Lyapunov exponents are computed with respect to a reference (or ?ducial) trajectory, and somehow generalize (see [Wiesel, 1999] for a discussion) the stability analysis performed in terms of the eigenvalues of the Jacobian for a ?xed point or of the logarithms of the eigenvalues of the transition matrix over one period for a periodic orbit (the so-called Floquet exponents). The standard working de?nition of a chaotic attractor [Parker & Chua, 1989] is that of an attracting bounded set with at least one positive exponent. Several methods have been proposed in the literature [Parker & Chua, 1989; Ramasubramanian & Sriram, 1999] to numerically compute the Lyapunov exponents, but the most interesting, both computational and theoretical, is based on the QR decomposition of the tangent map [Janaki et al., 1999;

713

Ramasubramanian & Sriram, 1999; Thi?eault, 2002]. In this paper we brie?y review this method and apply it to a general class of bilinear systems. Speci?c results are obtained for the buck converter [Deane & Hamill, 1990; Fossas & Olivar, 1996], although the same procedure should work for any other piecewise linear system. We write the differential equations satis?ed by the Lyapunov exponents and the ancillary angular variable, and solve them numerically, using the dominant attracting set as the ?ducial trajectory, for each value of the bifurcation parameter. We also write the equations when a periodic orbit is used as reference, and in this case we are able to solve them analytically, so that a closed-form expression is deduced for the highest Lyapunov exponent in this case, and the result coincides with the real part of the Floquet exponent for the same orbit [Batlle et al., 1999] (see also [Lim & Hamill, 1999] for a numerical computation of the Lyapunov exponents of the buck converter incorporating M¨ uller’s method [M¨ uller, 1995]).

714 C. Batlle et al.

Using adimensional variables (voltage in units of Vref , current in units of Vref /R and time in units of RC ) the di?erential equations of the system are (we use the same symbol t for the scaled time) x ˙ = ?x + y , (3)

y ˙ = ?γx + γV θ (r (t) ? x) ,

(a) The buck converter.

where v → x, i → y , V = E/Vref and γ = R2 C/L, and r (t) has been correspondingly scaled. In all the numerical results presented in this paper, the values Vref = 11.3 V, R = 22 ?, C = 47 · 10?6 C, L = 0.02 H, Vu = 8.2 V, Vl = 3.8 V, T = 400 · 10?6 s and a = 8.4 are used [Deane & Hamill, 1990]. The buck converter is a special case of a periodic bilinear system, a system switching between two linear vector ?elds, which for our purposes we take as x ˙ = Ax + a if α, x + h(t) > 0 , Bx + b if α, x + h(t) < 0 . (4)

(b) The periodic ramp governing the switch. Fig. 1. The basic scheme of the buck converter and its ramp.

Figure 1 represents the basic scheme of the buck converter. When it operates in continuous current mode (i.e. when the current through the inductor never drops to zero) it switches between two topologies according to whether v (t) is above or below a periodic ramp. This system exhibits a perioddoubling route to chaos using E as bifurcation parameter [Deane & Hamill, 1990]. Periodic orbits can be explicitly computed by solving the two linear systems and connecting the solutions at the crossing point, although a trascendental equation must be solved numerically [Fossas & Olivar, 1996]. We will see that this numerical information, i.e. the value of the switching time, is the only required piece of nonanalytical information in our analysis. The switching function can be modeled in terms of the step function θ as u([v ], t) = θ (r (t) ? v (t)) , (1)

Here A and B are matrices, a, b and α are vectors, x is the state variable, , denotes the scalar product and h(t) is a T -periodic function. In compact form this can be written as x ˙ = Bx + b + θ ( α, x + h(t))[(A ? B )x + (a ? b)]. (5) The paper is organized as follows. Section 2 presents the essential results about the QR decomposition and its use for the computation of the Lyapunov exponents. Section 3 treats the equations for a general two-dimensional bilinear system, solves them for the buck converter by numerical integration and computes an exact recurrence relation for T -periodic orbits, the limit of which is then computed analytically in Sec. 4. Section 5 deals with the extension to higher periodic orbits and Sec. 6 summarizes the results and points to some future lines of work.

2. The QR Decomposition

Let us consider a system given by a ?rst-order differential equation in Rn dx = F (x, t) , dt and let x0 (t) be a solution. x ∈ Rn , (6)

where [v ] denotes the functional dependence and r is the periodic ramp r (t) = Vref + with t computed mod T . Vl Vu ? V l + t, a aT (2)

Lyapunov Exponents for Bilinear Systems 715

Let us consider a nearby trajectory and the di?erence between them, z = x ? x0 (t), and its evolution: dz ?F = dt ?x

x0 (t)

z ≡ DF(t) · z

(7)

1, . . . , n, are the diagonal elements of R(t), r (t) is R(t) with the ith row divided by ?i (t), and d(t) is the lower-triangular matrix that a?ects the Gram– Schmidt orthogonalization of r (t), it can be shown [Thi?eault, 2001a] that the Lyapunov exponents are given by λi = lim

t→∞

The solution to (7) is z (t) = M (t)z (0), where

t

M (t) = Te

0

DF (s)ds

(8)

1 (log ?i (t) ? log dii (t)) . t

(14)

and T represents the time-ordered exponential or Peano–Baker series [Rugh, 1996]. Geometrically, M (t) is the tangent map of x0 (0) → x0 (t). The SVD (singular value decomposition) allows the factorization of any matrix, and in particular of M (t), as M (t) = U (t) · A(t) · V (t)

T

The dii are bounded, so they are irrelevant in the limit and hence the Lyapunov exponents can be computed as the asymptotic value of log(? i (t))/t. The equation satis?ed by the tangent map M is ˙ = DF · M M (15) with M (0) = I. Using the QR decomposition and operating with QT and R?1 one gets ˙ + RR ˙ ?1 = QT · DF · Q ≡ S . QT Q (16)

(9)

where both U (t) and V (t) are orthogonal (V T V = I = U T U ), and A(t) is diagonal. Since M T (t)M (t) = V (t)A(t)AT (t)V T (t) (10)

˙ ?1 is upper-triangular, with diagoNotice that RR ˙ i ??1 . The orthogonal manal elements given by ? i trix Q can be parameterized as [Janaki et al., 1999] Q = Q(1,2) Q(1,3) · · · Q(1,n) Q(2,3) · · · Q(n?1,n) , (17) where Q(i,j ) is a rotation of angle φij in the (i, j ) plane. Hence, M is parameterized by n(n ? 1)/2 angles, plus the n diagonal elements ? i of R, plus the strictly upper n(n ? 1)/2 elements of R, the latter playing no role in the computation of the Lyapunov exponents, as we are about to see. In˙ is antisymmetric, the diagonal deed, since QT Q elements in (16) yield ˙ i = Sii (φ)?i , ? i = 1, . . . , n . (18)

is semi-de?nite positive and symmetric, its eigenvalues are non-negative (the diagonal entries of A(t)AT (t)) and correspond to a basis of orthonormal eigenvectors (the columns of V (t)). The singular values of M (t) are the diagonal elements of A(t). The in?nite-time Lyapunov exponents associated with x0 (t), λi , are the logarithms of the eigenvalues of Λ = lim (M T (t)M (t)) 2t

t→+∞

1

(11)

whenever that limit exists. Although the SVD of the tangent map yields a nice splitting of the expanding and contracting spaces, it is the QR decomposition, also called continuous Gram–Schmidt orthogonalization [Ramasubramanian & Sriram, 1999], the one which provides a formulation with the minimal number of variables [Janaki et al., 1999]. In our case, it is given by M (t) = Q(t) · R(t) (12) where Q(t) is orthogonal and R(t) is uppertriangular and with positive diagonal elements. Notice that M (t)M (t) = R (t)R(t)

T T

The elements of S depend only on the angles θ and this is a separable equation, provided that the φij (t) are known. Comparison of the subdiagonal elements in (16) gives the equations ˙ )ij = Sij (φ) , (QT Q i>j. (19)

With the given parameterization of Q, the left-hand ˙ ij and the equation can be side of (19) is linear in φ rewritten as ˙ ij = gij (φ) , φ i >j. (20)

(13)

so it is R(t) which contains all the information about the Lyapunov exponents. If ?i (t), i =

This is a system of highly nonlinear equations in the n(n ? 1) angles, which must be solved and the result plugged in (18).

716 C. Batlle et al.

3. The QR Decomposition for a Bilinear System

For the vector ?eld given by (5) the Jacobian is ?F ?x ≡ DF(t) = B + θ ( α, x0 (t) + h(t))(A ? B ) + δ ( α, x0 (t) + h(t))[(A ? B )x0 (t) + (a ? b))αT ] , (21) where δ = θ is the Dirac delta function. The simplest case is when α, x0 (t) + h(t) has a single simple zero at t = tc mod T for each period. Assuming the (A, a) topology at the beginning of each period, we have θ ( α, x0 (t) + h(t)) = 1 ? θ (t ? tc ) δ ( α, x0 (t) + h(t)) = (22)

Computation of QT DFQ then yields (t is computed mod T unless otherwise speci?ed) ˙1 ? = f11 cos2 α ? (f12 + f21 ) sin α cos α ?1 + f22 sin2 α (29) ˙2 ? = f11 sin2 α + (f12 + f21 ) sin α cos α ?2 + f22 cos2 α (30) α ˙ = f12 sin2 α ? (f11 ? f22 ) sin α cos α ? f21 cos2 α , (31) where fij are obtained from the elements of A, B and C as fij = aij + (bij ? aij )θ (t ? tc ) + cij δ (t ? tc ). Notice that d log(?1 ?2 ) dt = f11 + f22 = a11 + a22 + (b11 + b22 ? a11 ? a22 )θ (t ? tc ) + (c11 + c22 )δ (t ? tc ) (32) and that this is the divergence of the vector ?eld given by the right-hand side of (5) computed on x0 (t). Using (t ? tc )δ (t ? tc ) ≡ 0, Eq. (32) can be integrated and we get log(?1 ?2 ) = (a11 + a22 )t + (b11 + b22 ? a11 ? a22 )(t ? tc )θ (t ? tc ) + (c11 + c22 )θ (t ? tc ) . (33) Since the equations for ?1 and ?2 are decoupled, this means that ?2 can be computed once ?1 is known. Equations (29) and (31) must be solved for each period [nT, (n + 1)T ], with initial conditions in each period given by the ?nal values in the preceeding period.1 Let us ?rst solve the equation for α. The presence of δ (t ? tc ) means that α(t) is not continuous at t = tc and that the corresponding jump must be computed; for t0 ≤ t < tc and tc < t ≤ t0 + T we can forget about the delta function. In the regime that we are considering, the eigenvalues of A and B must be complex with negative real part. Hence ?a ≡ (a11 ? a22 )2 + 4a12 a21 < 0 and a11 + a22 < 0 (34) and similarly for the elements bij of B .

x0 (t)

1 δ (t ? tc ). ˙ (tc )| | α, x ˙ 0 (tc ) + h (23)

Thus, de?ning the time-dependent matrix 1 D (t) = [(A ? B )x0 (t) ˙ (tc )| | α, x ˙ 0 (tc ) + h + (a ? b))αT ] , (24) we get the Jacobian associated with this kind of periodic orbit DF(t) = A + θ (t ? tc )(B ? A) + δ (t ? tc )C , (25) where the multiplying δ (t?tc ) allows us to introduce the constant matrix C = D (tc ) instead of D (t). The preceeding equations are valid for an ndimensional system. For a two-dimensional system we only have an angle φ12 (t) ≡ α(t) and two functions ?i (t), and thus Q(t) = cos α(t) sin α(t) ? sin α(t) cos α(t) 0 α ˙ ?α ˙ 0 ?1 (t) r12 (t) 0 ?2 (t) . , (26)

˙ = QT Q and

(27)

R(t) =

1

(28)

Notice that the going away or joining of nearby trajectories is caused only by the di?erent times at which they change the topology inside the period; the change at the end of a period a?ects all nearby trajectories at the same time.

Lyapunov Exponents for Bilinear Systems 717

? For t0 ≤ t < tc we have α ˙ = a12 sin2 α ? (a11 ? a22 ) sin α cos α ? a21 cos2 α (35) Notice that a12 = 0 under the assumed conditions. In terms of α(t0 ) ≡ α0 , the solution to this equation is tan α(t) √ √ a12 ??a ??a a11 ? a22 + · tan (t ? t0 ) = 2a12 2|a12 | 2|a12 | 2|a12 | a22 ? a11 + arctan √ tan α0 + 2a12 ??a . (36)

Since λ2 (t) can be obtained from λ1 (t), let us set λ(t) ≡ λ1 (t), with equation ˙ = f11 cos2 α ? (f12 + f21 ) cos α sin α + f22 sin2 α . λ (42) We must again consider three cases. ? t0 ≤ t < tc . Dividing the equations satis?ed by λ and α we get dλ a11 ? (a12 + a21 ) tan α + a22 tan2 α = . dα a12 tan2 α + (a22 ? a11 ) tan α ? a21 (43) The right-hand side of this equation can be converted to a rational integral and λ(α) can thus be obtained. The corresponding constant of integration must be determined from λ(α0 ) = λ(t0 ). ? At t = tc , the quotient of the terms proportional to δ (t ? tc ) in the equations for λ and α yields

dλ c11 cos2 α ? (c12 + c21 ) cos α sin α + c22 sin2 α = . dα c12 sin2 α ? c12 cos2 α ? (c11 ? c22 ) cos α sin α

? At t = tc we have a jump in the value of α. To compute it, we can disregard all the terms on the right-hand side of (31) except those proportional to δ (t ? tc ). We get α ˙ = (c12 sin2 α ? c21 cos2 α ? (c11 ? c22 ) cos α sin α)δ (t ? tc ) and integrating around tc

t+ c t? c

(44)

(37)

If we call H (α) the primitive of the right-hand side of (44) with respect to α, integration between + t? c and tc gives

? + ? λ(t+ c ) = λ(tc ) + H (α(tc )) ? H (α(tc )) .

c12 sin α ? c21

t+ c

2

cos2

dα α ? (c11 ? c22 ) cos α sin α

(45)

=

t? c

δ (t ? tc )dt = 1 .

(38)

If we denote by G(α) a primitive of the integrand of the left-hand side, we get

?1 ? α(t+ c ) = G (1 + G(α(tc ))) .

? For tc < t ≤ t0 + T we get the same equation than for t0 ≤ t < tc , replacing aij by bij . The initial + condition is now λ(α(t+ c )) = λ(tc ), and t must be replaced by t ? tc in the solution. We will specialize now the above equations for the adimensional buck converter. One has A=B= b= ?1 ?γ 0 γV 1 0 , , α= a= 1 0 0 0 , (46)

(39)

? For tc < t ≤ t0 + T we have the same equation as above, with bij instead of aij . The solution can be then obtained from (36) replacing t 0 with tc and α0 with α(t+ c ). Once α(t) is known, we can solve the equation for ?1 . Let us introduce the notation ?i (t) = e

λi (t)

.

(40)

and h(t) = ?r (t), and the Jacobian computed on (x0 (t), y0 (t)) is DF(t) = ?1 1 0 + δ (r (t) ? x0 (t)) ?γ 0 ?γV 0 0 .

Then, according to (14), the Lyapunov exponents are given by 1 λi = lim λi (t) . (41) t→∞ t

(47)

718 C. Batlle et al.

The explicit form of Eqs. (18) and (19) for this system is ˙1 ? = ? cos2 α + (γ ? 1) sin α cos α ?1 + γV δ (r (t) ? x0 (t)) sin α cos α ˙2 ? = ? sin2 α + (1 ? γ ) cos α sin α ?2 ? γV δ (r (t) ? x0 (t)) cos α sin α

2

2

1.5

(48)

1

0.5

(49)

0

α ˙ = ? sin α ? cos α sin α ? γ (1 + V δ (r (t) ? x0 (t))) cos 2 α

(50)

-0.5 1 1.5 2 2.5 3 3.5

˙ 1 ?2 +?1 ? ˙ 2 = ??1 ?2 and hence Notice that ? ? t ?1 (t)?2 (t) = e , since ?1 (0) = ?2 (0) = 1. In terms of λi (t) this is λ1 (t) + λ2 (t) = ?1 , t (51)

Fig. 2.

LLE of the dominant attracting set in terms of V .

so that the sum of the two Lyapunov exponents for the buck converter is equal to ?1. As mentioned in the general case, this is a direct consequence of the divergence of the vector ?eld of the adimensional buck converter being equal to ?1. We have numerically integrated the equations for λ(t), α(t) and the system variables (x 0 (t), y0 (t)) for a range of values of the bifurcation parameter V , using an arctan with a high coe?cient to approximate the step function and its derivative to approximate the delta function. The corresponding largest Lyapunov exponent (LLE) is displayed in Fig. 2, where it is seen that the LLE is ?0.5 in a wide region. As this is just a numerical integration, the system evolves to the dominant attracting set. Hence, a positive LLE corresponds here to the presence of a chaotic attractor. The zones of existence of such an attractor shown in Fig. 2 are in accordance with what is known about the buck converter [Fossas & Olivar, 1996]. For T -periodic orbits with a single switching time tc on each period, Eqs. (48) and (50) boil down to (29) and (31). In our case only c21 = ? γV ≡ ?ν |r ˙ (tc ) ? x ˙ 0 (tc )| (52)

α ˙ = sin2 α + sin α cos α + (γ + νδ (t ? tc )) cos2 α .

(54)

The denominator in (52) is the di?erence in slope between the ramp and the voltage at the crossing point, and this is the only information about the periodic orbit which enters the computations. As discussed for the general case, the only relevant issue for the integration of the above equations is the computation of the jump discontinuities at t = tc , i.e. Eqs. (39) and (45). One easily gets

? tan α(t+ c ) ? tan α(tc ) = ν .

(55)

and

? λ(t+ c ) ? λ(tc ) =

1 + tan2 α(t+ 1 c ) log ? . 2 2 1 + tan α(tc )

(56)

Using these results, the di?erential equations can be integrated over each period and a recurrence relation for the values αn = α(nT ) and λn = λ(nT ) can be computed. Putting rn = tan ?tc + arctan one gets the rational map rn+1 = A + B + rn , 1 ? AB ? Arn (58) γ ? 1/4, tan αn + 1/2 ? , (57)

is nonzero and the two relevant equations are ˙ = ? cos2 α + (γ + νδ (t ? tc ) ? 1) sin α cos α , λ (53)

where A = tan ?T , B = ν/? and ? = while for λn

Lyapunov Exponents for Bilinear Systems 719

λn+1 = where

(tan2 αn+1 + 1)(tan2 αn + tan αn + γ )((Mn + ν )2 + Mn + ν + γ ) 1 1 log ? T + λn , 2 + M + γ) 2 (tan2 αn + 1)(tan2 αn+1 + tan αn+1 + γ )(Mn 2 n Mn = ?rn ? 1/2 . (60)

(59)

Figure 3 shows the results of an iteration of (59) for a typical value of V in the T -periodic regime. It is seen that λn behaves linearly with a bounded oscillation on top of it, i.e. λ(t) = λ · t + φ(t) with φ(t) bounded. Hence λ(t) . (62) t Figure 4 shows the asymptotic values of α n and λn /(nT ) for a range of values of the bifurcation parameter. √ In fact, it will be shown that, for B ≥ B ? = 2(1 + 1 + A2 )/A, the map (58) has an attracting ?xed point and (αn ) →? ?π/2 (depends slightly on V ), while the sequence (α n ) does not converge for B < B ? . From the ?gure it is seen that this corresponds to λn /(nT ) going to ?1/2 in the case of no ?xed point, and to a value which depends on V when there is a ?xed point. The precise form of the dashed curve in Fig. 4 is in accordance with the exact computation of the Floquet exponents for the periodic orbits (see e.g. Fig. 3 in [Batlle et al., 1999]). λ = lim

t→+∞

the results that in the previous section have been obtained by numerical iteration. This brings to light some interesting features of the rational map (58), which we rewrite as f (x) = A+B+x . 1 ? AB ? Ax (63)

(61)

The solutions to f (x) = x are x=? B ± 2 B2 A + B ? , 4 A (64)

and ?xed points do exist if the condition B2 A + B ? ≥0 4 A holds i.e. B ≥ B ? = 2(1 + 1 + A2 )/A (66) (65)

as previously stated. Condition (66) implies that √ 2 + 2 1 + A2 ν≥ ? = 10.22531653 , A (67)

4. Analytical Results for the Buck Converter

For the buck converter, the explicit knowledge of the maps (58) and (59) allows to prove analytically

20

The ?rst value of V for which ν veri?es (67) is V = 2.13 (this corresponds to E = 24.1 V). This result is in accordance with the numerical results

0

-20

-40

-60

-80

-100

20

40

60

80

100

120

140

160

180

200

Fig. 3.

Evolution of λn with (59).

Fig. 4. Asymptotic value of angle (solid) and LLE (dashed) for the T -periodic orbit in terms of V .

720 C. Batlle et al.

that we can see in Fig. 4. Only one of the two values given by (64), the one given by the “+” sign, yields a stable ?xed point. Setting tan2 αn + 1 , Cn = tan2 αn + tan αn + 1 Dn = (Mn + ν )2 + Mn + ν + γ . 2+M +γ Mn n (68) (69)

with D? = (M? + ν )2 + M? + ν + γ , 2+M +γ M? ? (77)

we can rewrite (59) as λn+1 = and hence λm Cm 1 m 1 + log(Dm?1 · · · D0 ) + λ0 ? T . = log 2 C0 2 2 (71) Therefore, λm λ = lim m→+∞ mT 1 1 log Cm = ? + lim 2 m→+∞ 2mT + lim 1 log(Dm?1 · · · D0 ) . m→+∞ 2mT (73) (72) Cn+1 1 1 1 log + log Dn + λn ? T . (70) 2 Cn 2 2

M? = ?r? ? 1/2 and r? the attracting ?xed point of (63). Notice that this depends on the bifurcation parameter V through ν . Computing (76) for values of V above 2.13 yields values coincident with the increasing part of the dashed line in Fig. 4. When (αn ) is not convergent, i.e. V < 2.13 a (nontrivial) computation shows that 1 log(Dm · · · D1 ) = 0 m→+∞ 2mT lim (78)

and so the LLE is independent of the bifurcation parameter in this range and equal to ?1/2, as was previously seen in the numerical iteration.

5. Higher Periodic Orbits

In this section we want to show how the computations done for the T -periodic orbit can be generalized for higher periodic orbits. For a 2T -periodic orbit with switching times t1 and t2 (one on each period of the ramp) one has r (t1 ) ? x0 (t1 ) = 0 and r (t2 ) ? x0 (t2 ) = 0, 0 < t1 < T < t2 < 2T mod 2T and so δ (r (t) ? x0 (t)) = 1 δ (t ? t 1 ) |r ˙ (t1 ) ? x ˙ 0 (t1 )| + 1 δ (t ? t2 ) . (79) |r ˙ (t2 ) ? x ˙ 0 (t2 )|

Since log[(x2 + 1)/(x2 + x + 1)] is bounded for x ∈ R, the ?rst limit in (73) is trivially zero. Let us consider the second one 1 log(Dm?1 · · · D0 ) m→+∞ 2mT lim = 1 log(Dm · · · D1 ) . m→+∞ 2mT lim (74)

The system to solve in [0, 2T ] is α ˙ = sin2 α + sin α cos α + (γ + ν1 δ (t ? t1 ) ˙ = ? cos2 α + (γ + ν1 δ (t ? t1 ) λ + ν2 δ (t ? t2 ) ? 1) sin α cos α , (81) + ν2 δ (t ? t2 )) cos2 α (80)

Assume ?rst that (αn ) has a limit, i.e. the map (63) has a ?xed point. Then the sequences (r n ), (Mn ) and (Dn ) are convergent too. Let us call D? the limit of (Dn ). Using the fact that the limit of the arithmetic means of a convergent sequence equals the limit of the latter, one easily has in this case log D? 1 log(Dm · · · D1 ) = lim m→+∞ 2mT 2T and hence, for the case of ?xed point, 1 log D? , λ=? + 2 2T (76) (75)

where ν1 = γV /|r ˙ (t1 )?x ˙ 0 (t1 )|, and ν2 = γV /|r ˙ (t2 )? x ˙ 0 (t2 )|, and thus the information about the orbit only enters through t1 and t2 . The above system can be solved analytically and iterated, using the ?nal values of an iteration as the starting ones of the next, and a recurrence relation for the values αn = α(2nT ) and λn = λ(2nT ) can be computed.

Lyapunov Exponents for Bilinear Systems 721

2.5

2

1.5

1

0.5

0

-0.5

-1

1

1.5

2

2.5

3

3.5

Fig. 5. Asymptotic value of LLE (solid) of the dominant attracting set, LLE (dashed) for the T -periodic orbits and LLE (short dash) for the 2T -periodic orbits in terms of V .

Figure 5 shows the LLE obtained with a direct numerical integration of the system and the LLEs resulting from the recurrence relations for α n and λn for both the T -periodic orbit and the 2T -periodic orbit (genuine 2T -periodic orbits appear only when the T -periodic ones become unstable) for a range of values of the bifurcation parameter V . Whenever a periodic orbit is stable, the LLE obtained from the corresponding recurrence coincides with the one coming from the numerical integration. Notice the existence of a small chaotic window just before the T -periodic orbit becomes unstable. As in the T -periodic case, analytical expressions for the LLE can be obtained for the 2T periodic case, and the whole procedure can be replicated for kT -periodic orbits. In this way, using increasing values of k , one can obtain analytical expressions for successive pieces of the LLE curve in Fig. 2, as long as the dominant attractor is a periodic orbit, i.e. before the onset of chaos. However, since we have an analytical expression for the LLE associated with a periodic orbit, be it stable or not, our results could be a departure point for obtaining approximate analytical expressions of the LLE in chaotic regime, in the spirit of periodic orbit theory [Cvitanovi? c et al., 2000]. Our procedure can be, in principle, easily extended to any piecewise linear system and, in particular, to other PWM-controlled converters.

arise from the QR decomposition of the tangent map with respect to a reference trajectory. We have numerically solved the equations and computed the largest Lyapunov exponent for a range of values of the bifurcation parameter, and analytically when the reference trajectory is a periodic one, either stable or unstable. When a periodic orbit is the dominant attractor, the analytical result coincides with the value obtained by numerical integration. The results also agree with what was previously known in the literature about the real part of the Floquet exponents of the periodic orbits of the buck converter. These computations have been extended to higher periodic orbits and we hope to be able to give analytical approximations to the LLE in the fully chaotic regime, using tools from periodic orbit theory [Cvitanovi? c et al., 2000]. It would also be interesting to compute higher order corrections to the linear LLE, based on the Hessian instead of the Jacobian [Thi?eault, 2001].

Acknowledgments

We would like to thank Profs. E. Fossas and G. Olivar for useful discussions. C. Batlle was partially supported by CICYT PB98-0920 and by DPI2000-1509-C03-02.

References

Batlle, C., Fossas, E. & Olivar, G. [1999] “Stabilization of periodic orbits of the buck converter by time-delayed feedback,” Int. J. Circuits Theor. Appl. 27, 617–631. Batlle, C., Fossas, E. & Olivar, G. [2000] “From Floquet exponents to control of chaos in piecewise linear systems,” Proc. IEEE Int. Symp. Circuits and Systems (ISCAS2000) II, pp. 100–103. Cvitanovi? c, P., Artuso, R., Mainieri, R. & Vatay, G. [2000] Classical and Quantum Chaos (webbook), Niels Bohr Institute, Copenhagen — available from www.nbi.dk/ChaosBook/. Deane, J. & Hamill, D. [1990] “Analysis, simulation and experimental study of chaos in the buck converter,” in IEEE Power Electronics Specialists Conf. 1990, Vol. II, San Antonio, Texas, pp. 491–498. Filippov, A. [1988] Di?erential Equations with Discontinuous Righthand Sides (Kluwer Academic Pub., NY). Fossas, E. & Olivar, G. [1996] “Study of chaos in the buck converter,” IEEE Trans. Circuit Syst.-I 43, 13–25. Janaki, T. M., Rangarajan, G., Habib, S. & Ryne, R. D. [1999] “Computation of the Lyapunov spectrum for continuous-time dynamical systems and discrete maps,” Phys. Rev. E60, 6614–6626.

6. Conclusions

We have written the di?erential equations for the Lyapunov exponents of the buck converter which

722 C. Batlle et al.

Lim, Y. H. & Hamill, D. C. [1999] “Problems of computing Lyapunov exponents in power electronics,” Proc. IEEE Int. Symp. Circuits and Systems (ISCAS1999) V, pp. 297–301. M¨ uller, P. C. [1995] “Calculation of Lyapunov exponents for dynamic systems with discontinuities,” Chaos Solit. Fract. 5, 1671–1681. Parker, T. S. & Chua, L. O. [1989] Practical Numerical Algorithms for Chaotic Systems (Springer-Verlag, NY). Ramasubramanian, K. & Sriram, M. S. [1999] “A comparative study of Lyapunov spectra with di?erent algorithms,” preprint nlin.CD/9909029.

Rugh, W. J. [1996] Linear System Theory (Prentice Hall, Upper Saddle River, NJ). Thi?eault, J.-L. [2001] “Di?erential constraints in chaotic ?ows on curved manifolds,” preprint nlin CD/0105010. Thi?eault, J.-L. [2002] “Derivatives and constraints in chaotic ?ows: Asymptotic behaviour and a numerical method,” Physica D172, 139–161. Wiesel, W.-E. [1999] “Stability exponents, separation of variables, and Lyapunov transforms,” preprint nlin.CD/9905002.