Computation of the real logarithm for a discrete-time nonlinear system

Systems & Control Letters 59 (2010) 33–41

Contents lists available at ScienceDirect

Systems & Control Letters journal homepage: www.elsevier.com/locate/sysconle

Computation of the real logarithm for a discrete-time nonlinear system Laura Menini ∗ , Antonio Tornambè Dip. di Informatica, Sistemi e Produzione, Università di Roma Tor Vergata, via del Politecnico, 1-00133 Roma, Italy

article

info

Article history: Received 19 May 2009 Received in revised form 12 October 2009 Accepted 2 November 2009 Available online 24 November 2009

abstract For a given nonlinear discrete-time dynamical system, the problem is considered of finding, if any, a continuous-time nonlinear dynamical system such that the given system is its exact sampled-data representation. Constructive solutions to the problem are given by geometric tools such as symmetries. © 2009 Elsevier B.V. All rights reserved.

Keywords: Sampled-data systems Symmetries Poincaré–Dulac normal form

1. Introduction The problem of finding a discrete-time model of a nonlinear sampled-data system has been widely studied in the literature; see, for example, [1–5]. The derivation of the exact model is conceptually simple but computationally difficult; hence most of the existing results consider approximate discrete-time models and their relevance for control purposes. This paper considers the inverse problem: that of finding, if any, some continuous-time system whose exact sampling coincides with a given discrete-time system. For linear systems the problem has been dealt with recently in [6], taking into account multi-rate sampling schemes, whereas, for nonlinear systems, the problem has been considered in [7] from the mathematical point of view; the authors of [7] give conditions for the existence of a formal vector field, with possibly complex coefficients, describing the required continuous-time system. Here, on the contrary, the results are constructive, since two proposed algorithms allow one to compute the required continuous-time dynamical system, which is described by an analytic vector field with real coefficients. Apart from the conceptual interest of the problem considered here, we mention its practical relevance in identification (see the related works [8,9]): normally, system identification is carried out in the discrete-time framework, but it is known that the real system is a continuous-time one, whose description might be needed. As a by-product of the derivations carried out to solve this problem, a convenient method for the computation in closed form of the exact discrete-time model of a sampled-data system is given.

∗

Corresponding author. Tel.: +39 06 7259 7432; fax: +39 06 7259 7460. E-mail addresses: [email protected] (L. Menini), [email protected] (A. Tornambè). 0167-6911/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.sysconle.2009.11.003

For brevity, the exposition is referred to planar systems; Section 6 explains how to extend the results to higher-dimensional systems. The two algorithms reported in this paper are based on the knowledge of a symmetry of the given discrete-time system, having some special properties. Some concepts useful to apply the proposed algorithms, such as those of the first integral and inverse integrating factor for discrete-time systems, are introduced and compared with the corresponding concepts for continuous-time systems. Notations. σ (A) denotes the set of the eigenvalues of A, ◦ denotes function composition, ∂∂x applied to a scalar function stands for the row gradient, ∂∂x applied to a vector function stands ∂f

∂f

for the Jacobian matrix, and div(f ) stands for ∂ x1 + ∂ x2 . The 1 2 acronyms CT and DT stand for Continuous-Time and Discrete-Time, respectively. 2. Problem definition and results for linear systems Consider two vector functions f (x), F (x) ∈ R2 , and the two dynamical systems dx(t )

= f (x(t )),

x ∈ Uf ⊆ R2 ,

(1)

Dx(t ) = F (x(t )),

x ∈ UF ⊆ R2 ,

(2)

dt

T

where x = x1 x2 , Dh(t ) stands for h(t + δT ) with δT > 0 for any scalar or vector function h, and Uf and UF are open and connected subsets of R2 ; assume that Uf ∩ UF 6= ∅. For the sake of simplicity, U will denote either Uf or UF or Uf ∩ UF if systems (1) and (2) are to be jointly considered and, unless otherwise



34

L. Menini, A. Tornambè / Systems & Control Letters 59 (2010) 33–41

specified, it is assumed that all the functions are analytic in U. This implies that systems (1) and (2) have unique maximal solutions x(t ) = Φf (t , x0 ), t ∈ R+ , t sufficiently small to avoid finite escape times, and x(νδT ) = ΨF (ν, x0 ), ν ∈ Z+ , ν sufficiently small, respectively, from the initial condition x0 ∈ U at time t = 0; Φf and ΨF are the CT and DT flows associated with f and F , respectively. In the first part of the paper, it is also assumed that f (0) = 0 and F (0) = 0, unless otherwise specified, whereas in Section 5 the origin does not have a special role. System (2) is a sampled-data representation of system (1) if F (x) = Φf (δT , x), ∀x ∈ U∗ ⊆ U, for some positive sampling time δT ; in such a case ΨF (ν, x) = Φf (νδT , x), for all integers ν for which both sides are well defined. Known results relative to linear systems are contained in the following remark and illustrated in the example right after. Remark 1. If f = Ax and F = AD x, then Φf (t , x) = eAt x and ΨF (ν, x) = AνD x. System (2) is a sampled-data representation of system (1) if and only if there exists a real matrix A and δT > 0 such that AD = eAδT . From Theorem 1 of [10] (see, also, [11]), for a square real matrix AD , the equation AD = eAδT , with the requirement that A and δT are real, has a solution if and only if det AD 6= 0 and each Jordan block of AD that corresponds to an eigenvalue with negative real part occurs an even number of times; for AD of dimension 2, if AD has a negative eigenvalue λ, then AD = diag{λ, λ}. Clearly, in such a linear case, there is no loss of generality in taking δT = 1; considering the equation AD = eA , A can be termed as the ‘‘logarithm’’ of AD and AD as the ‘‘exponential’’ of A. The equation AD = eA can have (1) no real solution, (2) a unique real solution, (3) a countable set of real solutions, (4) an uncountable set of real solutions. Example 1. In this example, the general results given in [10] are illustrated (the techniques explained there can be usedhto compute i 1 −1

the matrix solutions for AD = eA ). (1) If either AD = AD =

h

0 1

i

1 0

−1 1

or

, no real solution of the matrix equation AD = eA

exists, in the first case because det AD = 0 (hence nor does a complex solution h exist) i and in the second case because the Jordan form of AD is (2) If AD =

h

2 1

−1 0

0 1

(in this case a complex solution exists).

i

1 2

, since the eigenvalues of AD are distinct positive

A real numbers (λ1 = 1, λ2 =  3), then AD = e has a unique

1

solution A =  21 2

ln 3 ln 3

1 2 1 2

ln 3

. (3) If AD = ln 3

h

0

−1

i

1 0

, since the

eigenvalues of AD are distinct complex conjugate numbers, with non-negative real part, " (ȷ, −ȷ), then AD = eA has a # countable set of solutions; in fact, A =

1

0

2

1

π (1 + 4ν)

− π (1 + 4ν) 2

AD = e for any integer ν . (4) If AD = A

is a solution of

0

h

1 0

i

0 1

, since the eigenvalues

of AD are identical positive real numbers (λ1 = λ2 = 1), then A A D = e has an uncountable  set of solutions; for instance, A = −2

 −2

ac + bd

ad − bc c 2 + d2 ad − bc

νπ

2

νπ

2

a 2 + b2

ad − bc ac + bd ad − bc

νπ

νπ

is solution of AD = eA for any integer ν

(including ν = 0) and for any reals a, b, c , d such that ad − bc 6= 0. Adopting the nomenclature used in the case of linear systems, it seems to be natural to give the following definition. Definition 1. If F (x) = Φf (δT , x), ∀x ∈ U, for some positive and real sampling time δT , then f is a logarithm of F and F is the exponential of f , with sampling time δT .

The aim of this paper is to give solutions to the following problem. Problem 1. (i) Given real f and δT > 0, find the real exponential of f , i.e., a real vector function F , such that F (x) = Φf (δT , x), ∀x ∈ U. (ii) Given a real F , find a real logarithm f of F , i.e., a real vector function f , such that F (x) = Φf (δT , x), ∀x ∈ U, for some real δT > 0. Part (i) of Problem 1 has a unique solution in the neighborhood of any point x at which f (x) is analytic, for sufficiently small δT ; to find the solution one needs to compute the flow of f , which can be difficult, in general. On the other hand, Part (ii) of Problem 1 does not need to have an (even complex) solution; when it has a solution, for fixed δT , then such a solution is not necessarily unique, even in the neighborhood of a point x at which F (x) is analytic. 3. Preliminaries Definition 2. Two scalar analytic functions h1 (x) and h2 (x) are



T

functionally independent if, letting h = h1 h2 , the Jacobian of h is nonsingular in some open set. Two vector functions f and g, analytic in some open set, and not identically zero, are colinear if f = α g, for some function α ;  equivalently, if det f g = 0 in some open set. From now on, the dependences on time τ ∈ R are omitted, if not necessary. Consider a vector function g (x) ∈ R2 and the corresponding CT dynamical system: dx dτ

= g (x),

x ∈ Ug ⊆ R2 ,

(3)

where Ug is an open and connected subset of R2 such that Uf ∩ UF ∩ Ug 6= ∅. Note that Ug does not necessarily contain the origin. In the rest of the paper U will denote Uf , Uf ∩ Ug , UF , UF ∩ Ug or Uf ∩ UF ∩ Ug depending on which system or systems are considered. By assuming g analytic in U, system (3) has a unique maximal solution x(τ ) = Φg (τ , x0 ), τ ∈ R+ , τ sufficiently small, from the initial condition x0 ∈ U at time τ = 0. The Lie derivative Lf h ∈ R of h by f is Lf h := ∂∂ hx f ; the CT-Lie bracket [f , g ] ∈ R2 of f and g is

[f , g ] :=

∂g ∂f f − g = Lf g − Lg f , ∂x ∂x

and the DT-Lie bracket bF , g c ∈ R2 of F and g is

bF , g c := g (F ) −

∂F g = g ◦ F − Lg F . ∂x

The definition of CT-Lie bracket is standard (see, e.g., [12–14]), whereas the DT-Lie bracket is used, without a formal definition of it, in [7]. The name DT-Lie bracket seems suitable because such an operation between two vectors has for DT systems the same meaning that the CT-Lie bracket has for CT ones [12, Section IV.7]. Note, however, that the well-known property [g , f ] = −[f , g ] of the CT-Lie bracket does not hold, in general, for the DT-Lie bracket: it holds if both F and g are linear, i.e., F (x) = AD x and g (x) = Cx for constant matrices AD and C ; in this case bAD x, Cxc = [AD x, Cx]. A first integral of the CT system (1) is a scalar function I (x) : U∗ → R, analytic in U∗ , such that Lf I = 0∀x ∈ U∗ , with U∗ being an open and connected subset of U; if I is a constant, then the first integral is said to be trivial; otherwise, it is non-trivial. A similar definition can be given for a first integral J (x) of system (3). On the other hand, for the DT system (2), one has to require I (F ) = I ◦ F = I. For brevity, a first integral of system (1) will be called a CT-first integral of f , and similarly for system (2).

L. Menini, A. Tornambè / Systems & Control Letters 59 (2010) 33–41

The flow y = Φg (τ , x) of system (3) qualifies as a oneparameter group of transformations; the orbits of the CT system (1) are mapped into themselves by y = Φg (τ , x), while preserving the time parameterization, if and only if

Φf (t , x) ◦ Φg (τ , x) = Φg (τ , x) ◦ Φf (t , x),

∀x ∈ U, ∀t , τ small enough.

(4)

Orbits of the DT system (2) are mapped into the same orbits by y = Φg (τ , x), while preserving the time parameterization, if and only if F (x) ◦ Φg (τ , x) = Φg (τ , x) ◦ F (x),

∀x ∈ U, ∀τ small enough. (5)

The following lemma is classical: for the CT case see Lemma 2.34 of [14]; the DT case can be proved in a wholly similar way. Lemma 1. The following relations hold: (4) holds ⇔ [f (x), g (x)] = 0, (5) holds ⇔ bF (x), g (x)c = 0,

∀x ∈ U, ∀x ∈ U.

g

∂J ∂x

T

`



−1

functions J1 and J2 such that J := J1 J2 satisfies = f g ; this implies that Lf J2 = Lg J1 = 0 and Lf J1 = Lg J2 = 1; hence J2 is a first integral of f and J1 is a first integral of g. Conversely, if there exist two functionally independent functions J1 and J2 such



that Lf J1 = 1 and Lf J2 = 0, letting J := J1 first column of


∂ J −1 ∂x


∂ J −1 ∂x

(6) `

the condition Ax, x11 x22 ej = 0, one obtains

are closed; i.e., there exist two functionally independent



`

where ej is thehjth column of i the 2 × 2 identity matrix; hence, from

Remark 2. In the CT case the notions of symmetry and first integral are strictly related. If g is a CT-symmetry of f (hence, f is a CT-symmetry of g), not colinear with f , then the two rows of matrix f

The aim of this subsection is to recall some background material about the Poincaré–Dulac normal form: for CT systems such a normal form is a classical tool (see, e.g., [16] or [17]), whereas for DT systems the reader is referred to [7]. Assume that f and F are analytic at x = 0, f (0) = 0 and F (0) = 0, and denote by Ax and AD x their linear parts, respectively; let f (x) = Ax + h(x) and F (x) = AD x + k(x), with h(x), k(x) having zero linear part. A square matrix is semi-simple if it can be rendered diagonal with a, possibly complex, linear change of coordinates. If A is semisimple, f is in the CT-Poincaré–Dulac normal form if [Ax, h] = 0. If AD is semi-simple, F is in the DT-Poincaré–Dulac normal form if bAD x, kc = 0. In the CT case, under the assumption that A is diagonal, A = diag{λ1 , λ2 }, one can write h(x) as the sum of a possibly infinite number of terms `

Definition 3. g is a CT-symmetry of f if [f , g ] = 0; g is a DTsymmetry of F if bF , g c = 0.

−1

3.1. The Poincaré–Dulac normal form for CT and DT systems

x11 x22 ej ,

The following definition follows from Lemma 1.



35

T

J2 , one has that the

coincides with f and the second column of

is a CT-symmetry of f .

The notions of symmetry and first integral play a crucial role for the operations of logarithm and exponential, since they are preserved under the operation of sampling; this statement is formalized in the following theorem, whose simple proof is omitted for brevity. Theorem 1. Under the assumption that F (x) = Φf (δT , x), for some positive real δT , then: (i) any CT-symmetry g of f is a DT-symmetry of F ; (ii) bF , f c = 0; (iii) any first integral I of (1) is a first integral of (2). Note that it is not necessarily true that, if bF , g c = 0 for some g, then g is a logarithm of F . As for the converses of Statements (i) and (iii), they hold under additional conditions, as stated in Lemma 3 reported in Section 4 and in its proof. For the next statements, let n ∈ Z, n ≥ 1. Given a vector of  T integers r = r1 . . . rn , define the integer dilation δεr x :=

T  r ε 1 x1 . . . ε rn xn , for any scalar real ε 6= 0. If ri = 1, i = 1, . . . , n, then δεr x is called the standard dilation. A scalar function h(x) : Rn → R is homogeneous of degree m ∈ Z, with respect to δεr x, if (see [15]) h(δεr x) = ε m h(x). A vector function f :=  T f1 . . . fn : Rn → Rn is homogeneous of degree m, with r respect to δε x, if fi is homogeneous of degree ri − m with respect to  T δεr x, i = 1, . . . , n. Let g = r1 x1 . . . rn xn ; if f is homogeneous of degree m with respect to δεr x, then [f , g] = mf .

  ` ` ` ` ` ` `1 λ1 x11 x22 + `2 λ2 x11 x22 ej = λj x11 x22 ej ; this last condition holds if and only if the following CT-resonance condition holds:

`1 λ1 + `2 λ2 = λj ,

(7)

for some non-negative integers `1 , `2 such that `1 + `2 ≥ 2; the corresponding monomial (6) is called a resonant term. In the DT case, under the assumption that AD is diagonal, AD = diag{µ1 , µ2 }, one can j write k(x)kas the sum of terms (6); hence, `

`

`

`

from the condition AD x, x11 x22 ej

= 0, one obtains

(µ1 x1 )`1 (µ2 x2 )`2 ej = µj x11 x22 ej ; this last condition holds if and only if the following DT-resonance condition holds: `

`

µ11 µ22 = µj ,

(8)

for some non-negative integers `1 , `2 such that `1 + `2 ≥ 2. Remark 3. Vectors h(x) and k(x) in the Poincaré–Dulac normal forms can be non-zero only if the eigenvalues of the linear part are such that the needed resonance condition holds. Note that the resonance conditions (7) and (8) are different; they are related in the following sense: if AD = eAδT , then, for each resonance condition satisfied by the eigenvalues of A, there is a corresponding resonance condition satisfied by the eigenvalues of AD ; the converse is not necessarily true since the eigenvalues of AD might satisfy, for some values of δT , more resonance conditions than those of A. Any f can be formally transformed into its CT-Poincaré–Dulac normal form through a formal series y = Q (x) to be iteratively computed by increasing step by step the order of Q (x); some conditions (see [16–18]) guarantee that such a formal series is ˜ , to an analytic diffeomorphism. If a converconvergent, in some U gence condition is not required, then by the Borel Lemma (see [19]), there exists a C ∞ -transformation such that the transformed f differs from its normal form for a flat vector function, i.e., by a vector function that is C ∞ in the neighborhood of the origin and has all the derivatives equal to 0 at x = 0; this also means that, for any arbitrarily high integer m > 0, there exists a polynomial diffeomorphism such that the transformed f differs from its normal form only for terms of order higher than m. Similar facts hold for the DT-Poincaré–Dulac normal form.

36

L. Menini, A. Tornambè / Systems & Control Letters 59 (2010) 33–41

3.2. Conditions on the linear parts of f and F

4. Computation of the real logarithm for ‘‘linearizable’’ systems

By considering a modified version of Problem 1, adding the requirement that both f and F are analytic at x = 0, with f (0) = 0, F (0) = 0, it is possible to derive the relation between the linear parts of F and of its logarithms. Actually the same relation could be derived under the weaker requirement that both f and F are smooth at x = 0, i.e., differentiable an arbitrary number of times. In the following, the requirement that functions f and F are analytic at x = 0 will be made either for simplicity, or when using results on normal forms. In what follows, when some function h isP analytic at x = 0, it +∞ i will be identified with the formal series h = i=0 Hi x , where the Hi ’s have suitable dimensions and

If there exists a change of coordinates y = Q (x) such that the system (2) in the new coordinates is linear, i.e.,

x0 = 1,

x1 =



x1 , x2

x31  x21 x2   x3 =  x1 x2  , . . . . 2 x32





x21 2 x = x1 x2  , x22

 



If f is analytic at x = 0, i.e., f (x) = A1 x + A2 x2 + A3 x3 + · · ·, with A1 = A, then, for small t, Φf (t , x) is also analytic at x = 0, i.e.,

Φf (t , x) = D1 (t )x + D2 (t )x + D3 (t )x + · · · , 2

3

with D1 (0) = E (E being the identity matrix), Di (0) = 0, i ≥ 2 and matrices Di (t ) analytic for small times. Then, from x˙ = f (x), one has

˙ 1 (t )x + · · · = A1 D1 (t )x + · · · , D which, taking into account the condition D1 (0) = I and the notation A1 = A, yields D1 (t ) = eAt . This yields the following lemma, valid for any dimension of the state space. Lemma 2. If f is analytic at x = 0, with f (0) = 0, then, for sufficiently small t ≥ 0, Φf (t , x) is also analytic at x = 0; if Ax is the linear part of f , then eAt x is the linear part of Φf (t , x). The following corollary, written for planar systems, gives a necessary condition for the solvability of Part (ii) of Problem 1. Corollary 1. If F is analytic at x = 0, with F (0) = 0 and with linear part AD x, a necessary condition for the existence of a real logarithm f of F , analytic at x = 0, is that det AD 6= 0 and, if AD has a negative eigenvalue λ, then AD = diag{λ, λ}. The discussion yielding Lemma 2 is based on the assumption that both F and f are analytic; it is important to stress that F being analytic at x = 0 does not imply that each of its logarithms f is analytic at x = 0, as shown in the following example. Example 2. Consider the CT system Sf described by f = q 4

1 x41

+ 2x22

  −x 2 x31

,

Dy(t ) = Q ◦ F ◦ Q −1 (y(t )) = F˜ (y(t )) = AD y(t ), for some matrix AD , then a logarithm of F is given by f (x) =



∂ Q (x) ∂x

 −1

AQ (x),

where A is a matrix, if any, such that AD = eAδT ; matrix A can be computed with the techniques in [10,11] recalled in Remark 1. Even when the existence of the diffeomorphism Q (x) that linearizes F can be guaranteed, its computation can be a difficult task. The algorithm described in this section can be used to avoid the computation of Q (x), and also to find a real logarithm for some F that cannot be linearized by a change of coordinates. Assumption 1. Let IF denote the set of DT-first integrals of F (x). If I1 and I2 belong to IF then they are not functionally independent. Assumption 1 is not very restrictive and implies that, given any non-trivial I ∈ IF , IF only contains all the functions of I. Lemma 3. Under Assumption 1, if F is the exponential of f , f analytic at x = 0, for some sampling time δT , then bF , g c = 0 ⇔ [f , g] = 0. Proof. In view of Theorem 1(i), only the implication bF , g c = 0 ⇒ [f , g] = 0 needs to be proved. In view of Assumption 1, f cannot be zero (otherwise F would be the identity and every function would be a DT-first integral of F ). In the neighborhood of a regular point of f (any neighborhood of x = 0 contains regular points of f ), f admits a non-trivial symmetry g¯ and non-trivial CTfirst integrals (see the discussion on page 2 of [21]). Let If denote the set of CT-first integrals of f ; the set of all symmetries of f can be parameterized as Sym(f ) = {C1 (I )f + C2 (I )¯g , I ∈ If }. On the other hand, any DT-symmetry g˜ of F can be written as g˜ (x) = K1 (x)f (x) + K2 (x)¯g (x), (because f and g¯ are not co-linear), from which g˜ (F ) = K1 (F )f (F )+ K2 (F )¯g (F ). In view of Theorem 1(i) and (ii), since f (F ) = ∂∂ xF f and g¯ (F ) = ∂∂ xF g¯ , one has g˜ (F ) = K1 (F )

∂F ∂F f + K2 ( F ) g¯ . ∂x ∂x

(9)

Since g˜ is a DT-symmetry of F , then g˜ (F ) = ∂∂ xF g˜ , whence g˜ (F ) = ∂F K f + K2 g¯ ), from which, K1 and K2 being scalar, ∂x ( 1 g˜ (F ) = K1 (x)

∂F ∂F f + K2 ( x ) g¯ . ∂x ∂x

(10)

Note that ∂∂ xF has rank 2 in the neighborhood of x = 0, because F is the exponential of an f analytic at x = 0; hence, by Lemma 2, ∂ F (0) AD = ∂ x satisfies AD = eAδT , with A being the matrix of the linear part of f . Since rank ∂∂ xF = 2, by comparing (9) and (10), one obtains K1 (x) = K1 (F (x)), and K2 (x) = K2 (F (x)), i.e., K1 and K2 are DT-first integrals of F . Therefore, comparing the parameterization of Sym(f ) with the expression of g˜ (x), one has that the set of DTsymmetries of F coincides with Sym(f ) if the set IF of the DT-first integrals of F coincides with If . Under the assumptions made, this is true because I ∈ If ⇒ I ∈ IF , by Theorem 1(iii), and both If and IF only contain all the functions of I, I being any non-trivial function in If . 




which is not analytic at x = 0, being not even defined there. A



T

symmetry g of f is g = x1 2x2 ; therefore, using the algorithm outlined subsequently in Section 5.2, the exponential F of f , for some sampling time, can be computed. In particular, using the quasi-polar coordinates introduced by Lyapunov in his thesis (available, e.g., in [20]), considering the weights (1, 2) of the symmetry g, it can be seen that all the orbits of Sf are closed; the trajectories are periodic, all with the same period, which is independent of the initial condition. This means that there exists a positive δT (the period of such trajectories) such that Φf (t , x) = Φf (t + δT , x), for any t and x. Then, the exponential of f , with sampling time δT , is F (x) = Φf (δT , x) = Φf (0, x) = x; this reasoning shows that there are linear DT systems that are the exponential of some nonlinear CT system, the last being not analytic nor even defined at x = 0.

Assumption 2. F (x) admits a DT-symmetry g (x) such that g (x) is ∂ g (0) analytic at x = 0, g (0) = 0, ∂ x = C , σ (C ) = {r1 , r2 }, with r1 , r2 positive integers. Remark 4. Two different facts show that Assumption 2 is not too restrictive. First, it can be proven that, if f (x) can be linearized by

L. Menini, A. Tornambè / Systems & Control Letters 59 (2010) 33–41

means of a change of coordinates analytic at x = 0, then it admits a symmetry g (x) as in Assumption 2 with r1 = r2 = 1; in fact, if the near-identity change of coordinates y = Q (x) is such that F˜ (y) = AD y, then, in the new coordinates, g˜ (y) = y is a DT-symmetry of F˜ ; hence in the given coordinates F satisfies Assumption 2 with its


−1

symmetry g (x) = ∂∂Qx Q (x), with r1 = r2 = 1. Secondly, without assuming the existence of a change of coordinates linearizing F , assume that f is a logarithm of F with linear part Ax, with A being semi-simple with real eigenvalues, and that gˆ is a DT-symmetry of F , with linear part Cˆ x, with Cˆ 6= γ A, ∀γ ∈ R, Cˆ being semi-simple with real eigenvalues; then, the real numbers c1 and c2 in the expression g = c1 f + c2 gˆ can be chosen so that Assumption 2 holds. The key point for proving this fact is that matrices A and Cˆ commute; hence they can be rendered diagonal with a common change of coordinates. Matrices A and Cˆ commute because, in view of Lemma 2, Cˆ = eAδT , and it is well known that A and eAδT commute; see [11]. The DT-symmetry g (x), satisfying the conditions in Assumption 2, can be transformed by means of a near-identity analytic diffeomorphism z = Qg (x) into its Poincaré–Dulac normal form. In the remainder of this section it is assumed that g (x) is already in its Poincaré–Dulac normal form. As a consequence, if r1 = r2 , g¯ (x) = (1/r1 )g (x) = x is also a DT-symmetry; in such a case, F (x) is homogeneous of degree 0 with respect to the standard dilation, and therefore F (x) is linear. Since this case can be dealt with directly, let r1 6= r2 . Lemma 4. Under Assumptions 1 and 2, if g (x) is in Poincaré–Dulac normal form, a necessary condition for the existence of a real logarithm f of F , f analytic at x = 0, is that bF (x), Cxc = 0. Proof. If F is the exponential of f , then, in view of Lemma 3, g is a CT-symmetry of f . Since g (x) is in Poincaré–Dulac normal form, then g (x) = Cx + k(x), [k(x), Cx] = 0. By the results in [17] (that apply to CT systems), since f is analytic at x = 0 and [f , g ] = 0, then [f (x), Cx] = 0. Hence, by Theorem 1, bF (x), Cxc = 0.  Lemma 5. Under Assumptions 1 and 2, if r1 6= r2 , a necessary condition for the existence of a real logarithm f of F , f analytic at x = 0, is that the dynamic matrix of the linear part of F is diagonal, with positive eigenvalues. The technical results above motivate the following algorithm. Algorithm 1. This algorithm describes the computation of a real logarithm f , analytic at x = 0, of F , analytic at x = 0, under Assumptions 1 and 2, with r1 6= r2 and g (x) in Poincaré–Dulac normal form. Step 0. If F is linear, then compute a logarithm f of F by means of the results recalled in Remark 1. Stop. Step 1. Check that bF (x), Cxc = 0, otherwise Stop (by Lemma 4, F is not the exponential of any f analytic in x = 0). Step 2. Find a linear change of coordinates x˜ = Px such that C˜ = PCP −1 is diagonal, C˜ = diag(r1 , r2 ), with r1 ≤ r2 . Compute F˜ (˜x) = PF (P −1 x˜ ). Note that, since F˜ is homogeneous of degree 0 with respect to δεr x, and is analytic at x = 0, necessarily r2 = ν r1 for some integer ν > 1 (otherwise F would have been linear, and this case is excluded in this step). Step 3. Let

f˜ (˜x) =

37





a1 x˜ 1 , a2 x˜ 2 + a3 x˜ ν1

where, if b2 6= bν1 , then a1 = a3 =

1

ln(b1 ),

δT

b3



b2 − bν1

a2 =

1

δT

ln(b2 ), (11)

 ν ln(b2 ) − ln(b1 ) , δT δT 1

whereas, if b2 = bν1 , then a1 =

1

ln(b1 ),

δT

ν ln(b1 ), δT

a2 = ha1 =

a3 =

b3

δT bν1

.

(12)

Note that, when the resonance condition b2 = bν1 holds, then, if F is not linear, a change of coordinates, analytic at x = 0, that linearizes F does not exist; however, the system can be rendered linear by means of state immersion, similarly to what can be done for CT systems (see [22]). The formulae given at Step 3 can be derived by computing the flow Φf˜ (δT , x˜ ), for example, by the techniques in [22]. Remark 5. It is useful to clarify that this section deals with some simple cases in which either it is possible to compute a diffeomorphism y = Q (x) that linearizes the given F , or it is possible to compute a diffeomorphism z = Qg (x) that brings a known DTsymmetry g (x) of F (x) into its Poincaré–Dulac normal form. The next section deals with cases where it is not possible to compute such changes of coordinates to simplify the problem. Example 3. Consider system (2) and its DT-symmetry g (x) with

 F (x) =



3 3 1/3 x1 − x2 + x2 + x1 − x2




x2 + x1 − x2 3




2 3


2

 ,



x + 5x2 3 g ( x) = 1 . 2x2 An analytic diffeomorphism bringing g (x) in Poincaré–Dulac normal form g˜ (z ) is z = Qg (x); in the new coordinates z the given system is described by F˜ (z ) with





x1 − x2 3 , x2

Qg (x) =

g˜ (z ) =





z1 , 2z2

F˜ (z ) =

1/3z1 . z2 + z1 2





Algorithm 1 can be applied in the coordinates z, and, since matrix C is already diagonal, the change of coordinates at Step 2 is not needed. By the formulae at Step 3, taking δT = 1, one obtains the logarithm f˜ (z ), which can be transformed in the original coordinates: f˜ (z ) = f (x) =

− ln (3) z1 9/4 ln (3) z1 2





,

1/4 ln (3) x1 − x2 3

−4 + 27x2 2 x1 − 27x2 5
2 9/4 ln (3) x1 − x2 3









.

Alternatively, one can consider another DT-symmetry of the given system:



x1 + x2 2

 2x2 −

27 8

x1 − x2 3


2

   ,

  ˜F (˜x) = b1 x˜ 1 , b2 x˜ 2 + b3 x˜ ν1

g¯ (x) = 

which is the general form of F˜ analytic at x = 0 and homogeneous of degree 0 with respect to δεr x. If either b1 or b2 is negative or equal to zero, then Stop (by Lemma 5, F is not the exponential of any f analytic in x = 0); otherwise, a real logarithm f˜ of F˜ is given by

which, having linear part equal to Ex (E being the identity matrix), shows that the given F (x) can be actually linearized by a change of coordinates y = Q (x). The change of coordinates and the vector Fˆ (y) representing the system in the new coordinates are

 x2 −

9 8

x1 − x2


3 2

38

L. Menini, A. Tornambè / Systems & Control Letters 59 (2010) 33–41 3

" Q (x) =

x1 − x2
2 9 x2 + x1 − x2 3 8

#

" ,

Fˆ (y) =

1

y1

# .

3 y2

Clearly, one can easily compute the logarithm fˆ (y) of Fˆ (y); if δT = 1 is assumed, going back to the original coordinates, the same f (x) as above is obtained. 5. Computation of the real logarithm in the general case The goal of this section is to describe an algorithm to solve Problem 1(ii) for systems to which it is not convenient to apply the results in Section 4. The algorithm is based on some technical results reported in Sections 5.1 and 5.2. 5.1. The inverse integrating factor The notion of inverse integrating factor for the CT system (1) is classical: a function ω, not identically equal to 0, is an inverse integrating factor (a CT-inverse  integrating factor associated with f ) if the one-form ω1 f2 −f1 is closed, in which case there exists

In view of the importance of being able to compute several DTinverse integrating factors associated with a given F , some results are collected in the following lemma and in the discussion below, concerning semi-invariants and their computation. Lemma 6. (i) If ω1 and ω2 are two DT-inverse integrating factors ω associated with F , then I = ω1 is a DT-first integral associated with F . 2 (ii) If ω and I are, respectively, a DT-inverse integrating factor and a DT-first integral associated with F , then ω ˆ = ωI is a DT-inverse integrating factor associated with F . (iii) If ω1 and ω2 are two DT-inverse integrating factors associated with F , then ω = a1 ω1 + a2 ω2 is a DT-inverse integrating factor associated with F for any pair of constant real numbers a1 , a2 . Definition 5. A function s(x) is a semi-invariant of the DT system (2) (a DT-semi-invariant associated with F ) with characteristic function λ(x) if

Ds(x) = s(F (x)) = λ(x)s(x), and s(x) and λ−1 (x) do not have common zeros/poles in U.

The following definition of a DT inverse integrating factor is motivated by the subsequent Theorem 2.

The name semi-invariant comes from the fact that, if s is a semiinvariant, the set Is := {x ∈ U : s(x) = 0} is invariant for system (2). It is clear that any DT-inverse integrating factor is
a semiinvariant, with characteristic function λ(x) = det ∂∂Fx . For this reason, we report hereafter a computational result for polynomial F and s. Assume that s(x) is a polynomial DT-semi-invariant which is the linear combination of some polynomial functions pi , i.e.,

Definition 4. A scalar function ω, not identically equal to 0, is a DT-inverse integrating factor associated with F if

s = c1 p1 + c2 p2 + · · · + cm pm ,

a first integral I (x) of system (1) such that ∂∂xI = ω2 , ∂∂xI = − ω1 , 1 2 locally in some open and connected subset of R2 . It is well known (see, e.g., [23]) that a function ω, not identically equal to 0, is an inverse integrating factor of system (1) if and only if f

f

Lf ω = div(f )ω.

∂F ω (F ) = det ∂x 



ω (x) .

(13)

Let I be a non-trivial CT-first integral associated with f and ω be the corresponding CT-inverse integrating factor; then, f can be recast as

∂I f = ωS (1) ∂x

T

where S (c ) :=

h



,

0 c

(14) −c 0

i

, for any scalar c. Conversely, if f 6= 0 and

(14) holds, then I is a non-trivial CT-first integral associated with f and ω is the corresponding CT-inverse integrating factor. Theorem 2. Let f 6= 0 and (14) hold. If F (x) = Φf (δT , x), for some δT > 0, then the function ω appearing in (14) is a DT-inverse integrating factor associated with F . Proof. By statement (ii) of Theorem 1, bF , f c = 0, whence f (F ) = ∂F f ; also, using (14) (with x = F ), it can be shown that ∂x

ω(F ) S (1)



T  T ∂ I ∂F ∂I = ω S ( 1 ) . ∂ x x =F ∂x ∂x

(15)

Now, if I is a CT-first integral associated with f , it is also a DTfirst associated with F , i.e., I (F (x)) = I (x), which implies integral ∂I ∂F = ∂∂xI . Hence, from (15), ∂ x x=F ∂ x

ω(F ) S (1)



∂ I ∂x

T

x =F

T ∂F ∂F ∂ I =ω S (1) ∂x ∂x ∂ x x=F    T ∂F ∂ I ω(x) S (1) = det , ∂x ∂ x x=F 

T 

where the last equality is proved by verifying that DS (1) DT = det(D)S (1), for any 2 × 2 matrix D. Since ∂∂xI 6= 0, then ω(F ) =


det ∂∂Fx ω(x).



Ds = λs,

(16a) cm 6= 0,

(16b)

with the ci ’s being suitable reals. Consider the following matrix: p1  Dp1

p2 Dp2

 Γ :=  

.. .

.. .

Dm−1 p1

Dm−1 p2

··· ··· .. . ···

pm Dpm 



.. .

. 

Dm−1 pm

Theorem 3. If (16) hold and det(Γ ) is not identically 0, then s is a factor of det(Γ ). Remark 6. Theorem 3 can be extended, with proper amendments, to the non-polynomial case. In its application, good candidates to be DT-semi-invariants are all the factors of det(Γ ), or of its minors, and not only those that are linear combinations of p1 , . . . , pm , since Γ could be a minor of another matrix obtained with an enlarged choice of the functions pi . The proof of the following lemma is omitted for brevity; it is based on computations similar to those used in the proof of Theorem 2. Lemma 7. (i) If ω and I are, respectively, a DT-inverse integrating factor and a non-trivial DT-first integral associated with F , then g := ω S (1)



∂I ∂x

T

,

(17)

is a DT-symmetry of F ; moreover, ω and I are, respectively, a CTinverse integrating factor and a CT-first integral associated with g. (ii) If g is a DT-symmetry of F and I is a non-trivial first integral associated with both g and F (Lg I = 0 and I (F ) = I), then g can be written as in (17) for some ω that is a DT-inverse integrating factor associated with F . 5.2. A change of coordinates for the computation of the CT-flow The algorithm presented in the next section is based on the computation of the flow, which can be based on a change of

L. Menini, A. Tornambè / Systems & Control Letters 59 (2010) 33–41

coordinates as described hereafter. If [f , g ] = 0 and f and g are not colinear, consider the diffeomorphism y = J (x), where J is given in Remark 2. Note that such a change of coordinates usually is not analytic at the equilibrium points of the given dynamical system, whereas it is analytic in some open and connected set; usually the diffeomorphism y = J (x) is useful for studying the system in some open set having the origin, or another equilibrium of interest, as an accumulation point (see Example 4 below). In the coordinates y, f˜ (y) = e1 and g˜ (y) = e2 , with ei being the ith column of the 2 × 2 identity matrix E. The general solution of y˙ = f˜ (y) is y(t ) = e1 t + y(0). Therefore, in the original coordinates the general solution can be written as x(t ) = J −1 (e1 t + J (x0 )), if the initial condition is x(0) = x0 . Hence, the flow of (1) is Φf (t , x) = J −1 (e1 t + J (x)). This yields the following expression of the exponential of f with sampling time δT : F (x) = J −1 (e1 δT + J (x)) ; this formula represents a solution to statement (i) of Problem 1. 5.3. Computation of the real logarithm The algorithm proposed hereafter to compute a real logarithm of a given F is based on a DT-symmetry g of F , a DT-inverse integrating factor ω ˆ and a non-trivial DT-first integral ˆI associated with F such that Lg ω ˆ = div(g )ωˆ

Lg ˆI = 1.

and

(18)

The algorithm relies on Lemma 7, which allows the construction of a family of ‘‘candidate’’ logarithms fc (x) (parameterized by the real parameter c) that have ω ˆ and ˆI, respectively, as associated CT-inverse integrating factor and CT-first integral; the change of coordinates described in Section 5.2 (with f replaced by fc ) is used to verify that, for a proper choice of the parameter c, fc is actually a logarithm of F . Algorithm 2. This algorithm describes the computation of a real logarithm f of F , given a DT-symmetry g of F . Step 1. Find a TD-inverse integrating factor ω ˆ and a non-trivial TDfirst integral ˆI associated with F such that Eqs. (18) hold. Step 2. Write the vector function

∂ ˆI ∂x

fc = ω ˆ S (c )

!T ,

(19)

with c being a real number to be fixed later. Step 3. Find a family of diffeomorphisms y = Jc (x) (parameterized by the constant c) such that, in some open and connected domain, ∂ Jc ∂x

=



fc

g

 −1

. This is possible since, in view of conditions

−1

(18), [fc , g ] = 0, whence the rows of fc g are closed (see Remark 2). Step 4. Rewrite the given system (2) in the coordinates y = Jc (x), i.e., compute F˜c = Jc ◦ F ◦ Jc−1 . If F˜c is of the form



F˜c = y +

  δT 0

,

Remark 8. The computations needed at Step 1 can be carried out using Theorem 3 (and Remark 6) to find DT-semi-invariants of F , among which one can choose those that are actually DT-inverse integrating factors, and use them, on the basis of Lemma 6 to write more DT-inverse integrating factors. To write all possible DTfirst integrals associated with F , one can take into account that, since the system is planar, either Assumption 1 holds (with the implication mentioned right after it) or there exist two functionally independent DT-first integrals I1 and I2 associated with F , in which second case, given any pair of functionally independent I1 , I2 ∈ IF , IF is constituted by all the functions of I1 and I2 . Now, two examples are proposed to illustrate the application of Algorithm 2. It is stressed that in both cases Algorithm 1 cannot be applied, because the linear part of F is such that the Poincaré–Dulac normal form does not yield relevant information: in Example 4 the linear part of F is the identity matrix, hence the system is trivially in Poincaré–Dulac normal form, but this gives no hint on how to find coordinates in which the system has a simpler expression; in Example 5 the linear part of F is not semi-simple.

 Example 4. Let F =

for some real c = c¯ (δT ), then f := fc¯ (δT ) is a logarithm of F with sampling time δT > 0. Remark 7. If δT is not fixed, and there exists a logarithm f of F , then Steps 3 and 4 of the algorithm above can be skipped, since a logarithm is obtained by either c = 1 or c = −1 (at most one of the two values actually gives a real logarithm). On the other hand, since matrix S (c ) is linear, in view of the form of fc , it is clear that the actual value of the sampling time, assumed positive throughout the paper, affects the logarithm to be found just by a scaling of the time variable. This fact, which is evident in the formulae (11) and (12) for Algorithm 1, is confirmed by the examples below.



x1

x2 exp

1 − x1

T

2 − x1

x3 2 1

 T . A DT-

(1 − x1 )

x2 . It can be verified that, for x2 > 0, ωˆ = x21 x2 and ˆI = −x21 + ln x2 are a DT-inverse integrating factor



symmetry of F is g = 0

and a DT-first integral, respectively, associated with F , satisfying the requirements at Step 1 (a wholly similar analysis can be carried out for x2 < 0). The vector fc defined at Step 2 is fc =

 T −cx21 −2cx31 x2 , whereas the family of diffeomorphisms y = Jc (x) mentioned at Step 3 (in the neighborhood of any x such that x1 6= 0 and x2 6= 0) is given by y1 = cx1 , y2 = −x21 + ln x2 ; 1 the corresponding inverses are x1 =

1 , y1 c

x2 = e

y2 + 21

y c2 1

. The

given DT system is described in the new coordinates by F˜c =  T 1 1 y1 − y2 , which is in the form (20) for c = − δT . This shows c  T 1 2 2 3 that f = x1 x1 x2 is a real logarithm of F , with sampling δT δT time δT . It is easy to check that the linear parts of f and F , A = 0 and AD = E, respectively, satisfy the condition eAδT = AD of Lemma 2. Fig. 1 shows some trajectories of the DT system described by F (x) and of the CT system described by f (x). Note that F1 (x1 ) < 0 when x1 > 1, and this is not captured by the computed logarithm, which is valid in the region x1 < 1, x2 > 0. On the other hand, it can be seen that the CT system found as a logarithm has a finite escape δ time at tf = x T when x1 (0) = x1,0 > 0; for the flow Φf (δT , x) to 1,0

be well defined it is necessary that δT < tf , whence that x1 < 1.

 Example 5. Let F =

x1 + x2

x2

1 − x2

1 − x2

T . A DT-symmetry of F is

T

g = x2 0 . It can be verified that ω ˆ = x32 and ˆI = x2 1 are a DT-inverse integrating factor and a DT-first integral, respectively, associated with F , satisfying the requirements at Step 1. The vector



(20)

39

1+x

T

fc defined at Step 2 is fc = cx2 + cx1 x2 cx22 , whereas the family of diffeomorphisms y = Jc (x) mentioned at Step 3 (in the neighborhood of any x such that x2 6= 0) is given by y1 = − cx1 ,



2

1+x1 . The given DT system is described in the new coordinates x2 T 1 1 Fc . This y1 y2 , which is in the form (20) for c δT

y2 =

by ˜ =





+

=

c

 2 T

shows that f = δ1 x2 + x1 x2 x2 is a real logarithm of F , with T sampling time δT . Also in this case the linear parts of f and F satisfy the condition eAδT = AD of Lemma 2.



40

L. Menini, A. Tornambè / Systems & Control Letters 59 (2010) 33–41 10

corresponding to I1 and I2 (correspondence to be defined similarly to (18)), then a candidate real logarithm fc of F is given by

equilibrium points

8

fc = ω ˆc

6 x2



∂ I1 ∂x

T

 ×

2

ln(4)

 -3

-2

-1

0 x1

1

2

3

Fig. 1. Example 4: trajectories of the DT system (circles) and of its logarithm (continuous lines).

6. Extensions to higher-dimensional systems There are no conceptual obstructions to extending both Algorithms 1 and 2 to the higher-dimensional case, x ∈ Rn , with n > 2; however, since some amendments of the definitions, assumptions and conditions would be needed, and there are many possible cases to be covered, to wholly explain the results for n > 2 would result in a much longer paper; for this reason, just a sketch of application of the two extended procedures is given here, in order to illustrate the generality of the method. The key for extending Algorithm 1 can be understood by looking at its Step 3, where the given DT system written in local coordinates x˜ is either linear, or characterized by a vector function F˜ (˜x) analytic at x˜ = 0 and homogeneous of degree 0 with respect to a dilation δεr with positive integer weights r1 , r2 , . . . , rn ; the corresponding f˜ (˜x) is taken in the same class, namely homogeneous of degree 0 with respect to the same dilation δεr . Then, the formulae analogous to (11) and (12) can be computed case by case for each n-plet (r1 , r2 , . . . , rn ), following [22]. For instance, if n = 3, r1 = 1, r2 = r3 = 3, in the case of real and positive eigenvalues for the linear part of F , under the proper conditions and assumptions, the given F (x) can be expressed in local coordinates x˜ as



(21)

the corresponding f˜ (x) is taken as





a1 x˜ 1 3 f˜ (˜x) = a2 x˜ 2 + a3 x˜ 1  . a4 x˜ 3 + a5 x˜ 31

(22)

In the case when resonance conditions do not hold, i.e., b31 − b2 6= 0

and b31 − b4 6= 0, a real logarithm f˜ (˜x) of F˜ (˜x) is given by (22) with

  ln(b1 )

δT

,

a2 =

ln(b2 )

δT

b3 ln

,

a3 =

  a4 =

ln(b4 )

δT

,

(24)



x1

  δT    ln(2) 1 ln(32) 3   f (x) =  x + x a 2 3 1 ;  δ 62 δT   T   ln(4) 1 21

a5

δT

(25)

x31

it is easy to verify that such a f (x) can be written as in (24) with x1 ln(4) 1 ˆ = c = 12 δ and I1 = x3 − 63 a5 x31 , I2 =  2 and ω T

1 x2 − 62 a3 x31


3 3

1 a3 x1 being, respectively, the two DT-first integrals of F x2 − 62 and the DT-inverse integrating factor associated with F . When n > 3, the form (24) of the candidate logarithm fc (x) can be extended using the wedge product, which extends the vector product for dimension greater than three.

7. Conclusions The main results of this paper are two algorithms for the computation of a logarithm of a given DT system. Algorithm 1 can be used for a class of systems that can be linearized either by a change of coordinates or by a state immersion. On the other hand, Algorithm 2 can be used for a wider class of systems, including systems for which the Poincaré–Dulac normal form is trivial, but nonlinear, or cannot be defined, as shown in two examples. The two algorithms can be extended to the higher-dimensional case along the guidelines of Section 6. References



b1 x˜ 1 3 F˜ (˜x) = b2 x˜ 2 + b3 x˜ 1  ; b4 x˜ 3 + b5 x˜ 31

a1 =

T

where × denotes the vector product and c is a real number. As an example, if F has form (21) with b1 = 4, b2 = 2, b4 = 1, then a logarithm f of F is given by (22) and (23), i.e.,

4

0 -4

∂ I2 ∂x

b5 ln

,

a5 =

δT

b31

b31

b4

− b4

b31 b2

δT b31 − b2


, (23)


.

If some resonance condition holds, different formulae can be found by the techniques in [22]. Algorithm 2 can be extended by properly changing the form (19) of the candidate real logarithm fc of F . In particular, when the dimension of x is n = 3, if I1 and I2 are two DT-first integrals of F and ω ˆ is the DT-inverse integrating factor associated with F and

[1] J.W. Grizzle, P.V. Kokotović, Feedback linearization of sampled-data systems, IEEE Transactions on Automatic Control 33 (9) (1988) 857–859. [2] D. Ne˘sić, A.R. Teel, P.V. Kokotović, Sufficient conditions for stabilization of sampled-data nonlinear systems via discrete-time approximations, Systems & Control Letters 38 (1999) 259–270. [3] D. Ne˘sić, A.R. Teel, A framework for stabilization of nonlinear sampled-data systems based on their approximate discrete-time models, IEEE Transactions on Automatic Control 49 (7) (2004) 1103–1122. [4] D.S. Laila, A. Astolfi, Construction of discrete-time models for port-controlled Hamiltonian systems with applications, Systems & Control Letters 55 (2006) 673–680. [5] L. Grüne, D. Ne˘sić, J. Pannek, K. Worthmann, Redesign techniques for nonlinear sampled-data systems, Automatisierungstechnik 56 (2008) 38–48. [6] F. Ding, L. Qiu, T. Chen, Reconstruction of continuous-time systems from their non-uniformly sampled discrete-time systems, Automatica 45 (2) (2009) 324–332. [7] T. Gramchev, S. Walcher, Normal forms of maps: Formal and algebraic aspects, Acta Applicandae Mathematicae 87 (2005) 123–146. [8] S. Borchers, P. Rumschinski, S. Bosio, R. Weismantel, R. Findeisen, Model invalidation and system identification of biochemical reaction networks, in: Proc. 16th IFAC Symp. Identification and System Parameter Estimation, S. Malo, France, 2009. [9] M. Larsson, D.S. Laila, Monitoring of inter-area oscillations under ambient conditions using subspace identification, in: Proc. IEEE Power and Energy Society General Meeting, Calgary, Canada, 2009. [10] W.J. Culver, On the existence and uniqueness of the real logarithm of a matrix, Proceedings of American Mathematical Society 17 (5) (1966) 1146–1151. [11] F.R. Gantmacher, Theory of Matrices, 1st ed., vols. 1–2, Chelsea Publishing Co., New York, 1959. [12] W.M. Boothby, An introduction to differentiable manifolds and Riemannian geometry, Academic Press Inc., 1975. [13] A. Isidori, Nonlinear Control Systems, Springer, London, 1995. [14] H. Nijmeijer, A.J. van der Schaft, Nonlinear dynamical control systems, Springer-Verlag, 1990.

L. Menini, A. Tornambè / Systems & Control Letters 59 (2010) 33–41 [15] H. Hermes, Nilpotent and high-order approximations of vector field systems, SIAM Review 33 (2) (1991) 238–264. [16] J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer, 1983. [17] G. Cicogna, G. Gaeta, Symmetry and perturbation theory in nonlinear dynamics, in: Lecture Notes in Physics Monographs, vol. 57, Springer, 1999. [18] G. Gaeta, Poincaré normal and renormalized forms, Acta Applicandae Mathematicae: An International Survey Journal on Applying Mathematics and Mathematical Applications 70 (1) (2002) 113–131. [19] P. Hartman, Ordinary differential equations, in: Classics in Applied Mathematics, vol. 18, Society for Industrial and Applied Mathematics, 2002.

41

[20] A.M. Lyapunov, The general problem of the stability of motion, International Journal of Control 55 (3) (1992) 531–773 (special issue). [21] L. Menini, A. Tornambè, On the use of semi-invariants for the stability analysis of planar systems, in: Conference on Decision and Control, Cancun, Mexico, 2008. [22] L. Menini, A. Tornambè, Linearization through state immersion of nonlinear systems admitting Lie symmetries, Automatica 45 (8) (2009) 1873– 1878. [23] S. Walcher, Plane polynomial vector fields with prescribed invariant curves, Proceedings of the Royal Society of Edinburgh: Section A Mathematics 130 (2000) 633–649.