Chaotic dynamic and control for micro-electro-mechanical systems of massive storage with harmonic base excitation

Chaotic dynamic and control for micro-electro-mechanical systems of massive storage with harmonic base excitation

Chaos, Solitons and Fractals 39 (2009) 1356–1370 www.elsevier.com/locate/chaos Chaotic dynamic and control for micro-electro-mechanical systems of ma...

1MB Sizes 2 Downloads 50 Views

Chaos, Solitons and Fractals 39 (2009) 1356–1370 www.elsevier.com/locate/chaos

Chaotic dynamic and control for micro-electro-mechanical systems of massive storage with harmonic base excitation Manuel F. Pe´rez Polo a

a,*

, Manuel Pe´rez Molina

b,1

, Javier Gil Chica

a,2

Departamento de Fı´sica, Ingenierı´a de Sistemas y Teorı´a de la Sen˜al, Universidad de Alicante, Escuela Polite´cnica Superior, Campus de San Vicente, 03071 Alicante, Spain b Facultad de Ciencias Matema´ticas, Universidad Nacional de Educacio´n a Distancia. UNED, C/Boyero 12-1A, Alicante 03007, Spain Accepted 29 May 2007

Abstract This paper explores chaotic behaviour and control of micro-electro-mechanical systems (MEMS), which consist of thousands of small read/write probe tips that access gigabytes of data stored in a non-volatile magnetic surface. The model of the system is formed by two masses connected by a nonlinear spring and a viscous damping. The paper shows that, by means of an adequate feedback law, the masses can behave as two coupled Duffing’s oscillators, which may reach chaotic behaviour when harmonic forces are applied. The chaotic motion is destroyed by applying the following control strategies: (i) static output feedback control law with constant forces and (ii) geometric nonlinear control. The aim is to drive the masses to a set point even with harmonic base excitation, by using chaotic dynamics and nonlinear control. The paper shows that it is possible to obtain a positioning time around a few ms with sub-nanometre accuracy, velocities, accelerations and forces, as it appears in the design of present MEMS devices. Numerical simulations are used to verify the mathematical discussions.  2007 Elsevier Ltd. All rights reserved.

1. Introduction The relevance of nanotechnology is well recognized, so new developments and applications based on nonlinear dynamics, chaos, fractals, and nonlinear control, are nowadays reached in an interdisciplinary framework [1,2]. The term MEMS refers to mechanical microstructures (on the order to 10–1000 lm), such as sensors, valves, gears, gyroscopes and electric microactuators, fabricated using the same photolithographic processes used in the fabrication of CMOS devices [3–5]. The MEMS of massive storage consist of a large array of micrometric read/write heads called

*

Corresponding author. Tel.: +34 9 659 09673; fax: +34 9 659 09750. E-mail addresses: manolo@dfists.ua.es (Manuel F. Pe´rez Polo), [email protected] (M.P. Molina), gil@dfists.ua.es (J.G. Chica). 1 Tel.: +34 9 651 01988. 2 Tel.: +34 9 659 09673; fax: +34 9 659 09750. 0960-0779/$ - see front matter  2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2007.06.010

Manuel F. Pe´rez Polo et al. / Chaos, Solitons and Fractals 39 (2009) 1356–1370

1357

probe tips, which are held over/under a rectangular magnetic media sled where the data are stored and retrieved. The media sled has a linear motion, so by positioning the sled it is possible to transfer data to the probe tips. Details can be found in references [5–8]. The aim of this paper is to show that, from a nonlinear model, it is possible to obtain chaotic behaviour of a MEMS-based storage device, as well as analysing and designing efficient control laws for the positioning of the media sled. The mathematical model presented in this paper is formed by two masses connected by a hard spring and a linear damping to obtain two coupled Duffing’s oscillators [9–11]. The masses are upon a support-based which may be affected by external harmonic disturbance. The parameter values of the system are obtained from data reported in the bibliography [6]. From the analysis of one Duffing oscillator, and considering the Melnikov method, sufficient conditions for chaotic behaviour are obtained [12–16]. When the oscillators are completely actuated, it is researched how to design nonlinear control laws to drive the masses to an admissible equilibrium point. To do this end, the control law which give chaotic motion is commuted to a static feedback control law with constant force. The paper shows that it is possible to drive the media sled to the desired set point with minimum settling time, speed, acceleration and control forces [6,17]. The simulation values agree with the results published on experimental devices [4,17–19]. The previous control strategy is designed from a physical point of view, so a new systematic static feedback control law is researched from nonlinear control-based on differential geometry [20,21]. Finally, a discussion of the positioning and control laws based on very small time simulation is presented.

2. Statement of the problem In this paper a device of moving media with movement in only one dimension is considered (Fig. 1a). Description and details the media sled of these devices can be found in [6,7,17–19]. The mechanical model is formed by two masses connected by a nonlinear hard spring and a damping with the possibility of harmonic base disturbances. A schematic representation is shown in Fig. 1b. The control forces over the masses m1 and m2 are F1(t) and F2(t) respectively. L is the length of the spring without deformation, and the restoring force is defined by the following equation: F s ¼ K½ðr2  r1  LÞ þ as ðr2  r1  LÞ3 

ð1Þ

where K is the spring constant and as is a constant which defines a hard spring. From Fig. 1, the kinetic T, potential V and Rayleigh’s dissipation F functions are the following: 9     > = T ¼ 12 m1 r_ 21 þ r21 h_ 2 þ 12 m2 r_ 22 þ r22 h_ 2 h i ð2Þ   ; V ¼ K 12 ðr2  r1 Þ2 þ as ðr2  r1 Þ4 ; F R ¼ 12 B r_ 22  r_ 21 > 4

To simplify the mathematical model of the system the following auxiliary variables are introduced:  p1 ðtÞ ¼ r1 ðtÞ; p2 ðtÞ ¼ r2 ðtÞ  r1 ðtÞ  L ) p2 ðtÞ ¼ r2 ðtÞ  p1 ðtÞ  L p_ 1 ðtÞ ¼ r_ 1 ðtÞ; €p1 ðtÞ ¼ €r1 ðtÞ; €p2 ðtÞ ¼ €r2 ðtÞ  €r1 ðtÞ ) €r2 ðtÞ ¼ € p1 ðtÞ þ € p2 ðtÞ

ð3Þ

From the Lagrange equations, and considering Eqs. (2) and (3), the Lagrangian and the corresponding dynamical equations of the coupled masses can be written as follows [22]:    1 1  2 1 2 as 4 2 2 _2 2 _2 ð4Þ L ¼ m1 p_ 1 þ p1 h þ m2 ½ðp_ 1 þ p_ 2 Þ þ ðp2 þ p1 þ LÞ h   K p2 þ p2 4 2 2 2 ) m1 €p1 ðtÞ þ m2 ½€p1 ðtÞ þ €p2 ðtÞ  m1 p1 ðtÞh_ 2  m2 ½p1 ðtÞ þ p2 ðtÞ þ Lh_ 2 ¼ F 1 ðtÞ ð5Þ m2 ½€p1 ðtÞ þ €p2 ðtÞ  m2 ½p ðtÞ þ p ðtÞ þ Lh_ 2 þ K½p ðtÞ þ as p3 ðtÞ ¼ F 2 ðtÞ  Bp_ 2 ðtÞ 1

2

2

2

By removing the variables €p2 ðtÞ and €p1 ðtÞ of the first and second Eq. (5) respectively, and introducing the state variables: x1 ðtÞ ¼ p1 ðtÞ;

x2 ðtÞ ¼ p_ 1 ðtÞ;

x3 ðtÞ ¼ p2 ðtÞ;

x4 ðtÞ ¼ p_ 2 ðtÞ

the mathematical model of the system can be rewritten as:

ð6Þ

Manuel F. Pe´rez Polo et al. / Chaos, Solitons and Fractals 39 (2009) 1356–1370

1358

L2

Magnetic media

L1

Y

θ(t) = Aθ.cos(ωθt)

O

x1(t)

x2(t) – x1(t)

X

x2(t) = r2.cos(ωθt) Fig. 1. One-dimensional mechanical model of a MEMS device with harmonic base disturbances. Parameter values: m1 = 50 · 104 kg, m2 = 1 · 104 kg, K = 2156 N/m, as = 4.5 · 107 m1, B = 0.5 kg/s, L = 10 · 106 m, d = 9.486 (default value).

x_ 1 ðtÞ ¼ x2 ðtÞ x_ 2 ðtÞ ¼ mK1 x3 ðtÞ þ as mK1 x33 ðtÞ þ mB1 x4 ðtÞ þ x1 ðtÞh_ 2 þ um1 ðtÞ 1

9 > > > > =

> > > u1 ðtÞ u2 ðtÞ > K K 3 B 2 _ x_ 4 ðtÞ ¼  me x3 ðtÞ  as me x3 ðtÞ  me x4 ðtÞ þ ½x3 ðtÞ þ Lh  m1 þ me ; x_ 3 ðtÞ ¼ x4 ðtÞ

ð7Þ

where me = m1m2/(m1 + m2) is the reduced mass, and as we will see later, to simplify the mathematical treatment, the control forces have been modified as follows: ) ) F 1 ðtÞ 2 ðtÞ F 1 ðtÞ ¼ u1 ðtÞ þ mm11um  Fm2 ðtÞ ¼ um1 ðtÞ m1 e 1 1 ) ð8Þ 2 ðtÞ  Fm1 ðtÞ þ Fm2 ðtÞ ¼  um1 ðtÞ þ um2 ðtÞ F 2 ðtÞ ¼ mm11um e e e 1 1 Eq. (7) are two coupled-forced Duffing‘s oscillators, with coupling constants given by K/m1, asK/m1, B/m1 and K/me, asK/me, B/me respectively. Note that there is one-way coupling from the second oscillator to the first one, since the state variables x3(t) and x4(t) influence on to the velocity x_ 1 ðtÞand acceleration x_ 2 ðtÞ of the first oscillator, whereas the state variables x1(t) and x2(t) have no influence on the velocity x_ 3 ðtÞ and acceleration x_ 4 ðtÞ of the second one.

Manuel F. Pe´rez Polo et al. / Chaos, Solitons and Fractals 39 (2009) 1356–1370

1359

The study of chaotic motion of the coupled oscillators defined by Eq. (7) is carried out from the control signals u1(t) and u2(t), which may be used to decouple them. Let us suppose that there is no harmonic base excitation, i.e. h_ ¼ 0 and u2(t) = 0, and the following control law is applied: u1 ðtÞ ¼ K p x3 ðtÞ  uðtÞ Assuming that u(t) = A cos xdt, the equations of the second oscillator can be written as follows:  x_ 3 ðtÞ ¼ x4 ðtÞ x_ 4 ðtÞ ¼ ax3 ðtÞ  as x20 x33 ðtÞ þ e½c cos xd t  b0 x4 ðtÞ

ð9Þ

ð10Þ

where a ¼ x21  x20 , x21 ¼ K p =m1 , x20 ¼ K=me , b = B/me is a damping coefficient, eb 0 = b and ec = A/m1 with e  1. Now, taking Kp > K, we have that x21 > x20 and a > 0. Consequently, Eq. (10) are a Duffing’s oscillator in the form used to study chaotic behaviour [9–11,19]. By using Eq. (10) and from the Melnikov method [12–16], it is possible to deduce the necessary condition for chaotic motion, i.e: sffiffiffiffiffiffiffiffiffiffiffiffi   c 4 a3 pxd p ffiffi ffi cosh > f ðxd Þ; f ðxd Þ ¼ ð11Þ 3pxd 2as x20 2 a b0 This equation will be used in the analysis of the chaotic oscillation of the coupled oscillators in combination with the non-linear control laws.

3. Chaotic motion of the coupled oscillators In this section it is analyzed the dynamics of the system when a feedback law is applied to obtain chaotic oscillations. As we will see later, our purpose is to use the chaotic oscillations, in combination with nonlinear control laws, to improve the positioning of the system in a desired set point. To do this end, one could substitute the control law (9) into Eq. (7) to obtain: 9 x_ 1 ðtÞ ¼ x2 ðtÞ > >   > > K > = x_ 2 ðtÞ ¼ mK1  m1p x3 ðtÞ þ as mK1 x33 ðtÞ þ mB1 x4 ðtÞ  uðtÞ m1 ð12Þ > x_ 3 ðtÞ ¼ x4 ðtÞ > >   > > K ; x_ 4 ðtÞ ¼ m1p  mKe x3 ðtÞ  as mKe x33 ðtÞ  mBe x4 ðtÞ þ uðtÞ þ um2 ðtÞ m1 e where it is assumed that there is no harmonic base excitation i.e. h_ ¼ 0. It is not difficult to show that Eq. (12) cannot work properly. In fact, note that the equilibrium points of system (12) are deduced from the following equations: 9   K K > =  m1p x3e þ as mK1 x33e  mue1 ¼ 0 m1   ð13Þ Kp ;  mKe x3e  as mKe x33e þ mue1 þ um2ee ¼ 0 > m1 where the subindex ‘‘e’’ means equilibrium. However, Eq. (13) are cubic, so if one chose one value of Kp to obtain a real value of x3e from the first Eq. (13), it will be necessary to take the appropriate value of u2e to obtain the same equilibrium point from the second Eq. (13). Consequently, the control law (9) may be very cumbersome. To overcome this problem, it is possible to choose new functions of the signal control u1(t) and u2(t) defined as: ) u1 ðtÞ ¼ as K q p32 ðtÞ  B1 p_ 2 ðtÞ þ f1 ðtÞ ð14Þ u2 ðtÞ ¼ K p p2 ðtÞ  B2 p_ 2 ðtÞ þ f2 ðtÞ where the parameters Kp, Kq, B1 and B2 are undetermined and f1(t), f2(t) are new signals which may be defined to obtain chaotic behaviour. Substituting Eq. (14) into Eq. (7) and taking into account the discussion of the previous section to obtain a Duffing‘s oscillator, the following relations must be verified:   K Kp K me ¼  ) Kp ¼ K þ1 ð15Þ m1 m1 me me   Kq K K Kq K m1  ¼  ) Kq ¼ þ1 ð16Þ m1 m1 me m1 2 me

Manuel F. Pe´rez Polo et al. / Chaos, Solitons and Fractals 39 (2009) 1356–1370

1360

Eqs. (15) and (16) determine the values of Kp and Kq in Eq. (14). Parameters B1 and B2 may be determined as follows: Equalling the coefficients of variable x4(t) between second and four Eqs. (12)–(14) it is deduced that: ( B1 ¼ B þ m1 d B1 B B B1 B2   ð17Þ  ¼  þ ¼d) B2 ¼ 2me d  B 1  mm1e m1 m1 me m1 m1 where d is an equivalent damping coefficient whose value is chosen so that B1 is approximately equal to B. Note that m1 is very small and d  1. As we will see later, a wide range of values such as 5 6 d 6 2000 may be adequate. The signals of f1(t) and f2(t) are chosen to obtain chaotic behaviour i.e. f1 ðtÞ ¼ A1 cos xd 1 t; f2 ðtÞ ¼ A2 cos xd 2 t. From A1 = 0.05 N, A2 = (2me/m1)A1, xd1 = xd2 = 200 rad/s the values c/b 0 = 1.0542; f(xd) = 0.0215 are obtained. Consequently, the inequality (11) is verified, and the masses may reach chaotic oscillations. Since the system defined by Eqs. (7)–(14) is not autonomous, a new state variable has been added to avoid numerical problems, i.e: xd t ¼ x5 ðtÞ ) x_ 5 ðtÞ ¼ xd where the angular variable x5(t) is normalized between 0 and 2p. Fig. 2a shows a strange attractor in the phase plane defined by position and velocity of mass m2. The range of reachable points is between ±100 lm, which coincides with admissible values used in MEMS of massive storage. It is well known that it is very difficult to determine rigorously when a system has chaotic behaviour [12–15]. Therefore, indicators of chaotic dynamics such as Lyapunov exponents and power spectral density have been calculated, as shown in Fig. 2b and c. The power spectral density has a broad band with continuous energy vanishing, which is a typical feature of many chaotic systems [15]. The Lyapunov exponents have been calculated by linearizing Eqs. (7)–(14) throughout the trajectory and applying the Gram–Schmidt orthogonalization method [14,23,24]. The results show that there is a positive exponent which is an indicator of chaotic dynamics. Fig. 2d shows the time dependence of the Lyapunov exponent sum. This sum is very close to the divergence of the vector field defined by Eqs. (7)–(14) i.e: i¼5 X i¼1

NlðiÞ ¼ div ~ f ; div ~ f ¼

  i¼5 X ofi B B1 B2 ¼ d ¼  þ oxi me m1 m1 i¼1

ð18Þ

where the value d = 9.486 has been chosen, so the simulation process can be considered to be reasonably correct.

Fig. 2. Chaos indicators for the control law 14. Parameter values are indicated in legend of Fig. 1. Simulation values: A1 = 0.05 N, xd = 200 rad/s, integration time 1 s, simulation step 104 s. Parameter values of inequality (11): c/b 0 = 1.0542; f(xd) = 0.0215. (a) Strange attractor in the phase plane of mass m2. (b) Power spectral density. (c) Lyapunov exponents. (d) Sum of the Lyapunov exponents vs. time.

Manuel F. Pe´rez Polo et al. / Chaos, Solitons and Fractals 39 (2009) 1356–1370

1361

4. Positioning of masses with static nonlinear feedback laws In this section, the previous chaotic behaviour is used in combination with new nonlinear control laws. The purpose of these laws is to position the masses in an arbitrary set point inside of a reachable region filled by the strange attractor (see Fig. 2a). The problem of control with harmonic base excitations will be considered in the next section. The starting point will be Eq. (7) and the control signals u1(t), u2(t) defined by the equations: ) u1 ðtÞ ¼ as K q x33 ðtÞ  B1 x4 ðtÞ þ A1 cos xd t ð19Þ e u1 ðtÞ ¼ K p x3 ðtÞ  B2 x4 ðtÞ þ 2m A1 cos xd t m1 where the constants Kp, Kq, B1 and B2 are defined by Eqs. (15)–(17) respectively. Taking A1 = 0, the equilibrium point of system (7) is given by x1e, x2e = 0, x3e, x4e = 0, u1e ¼ as K q x33e , u2e = Kpx3e, being x1e and x3e a desired set point. From the values of x1e and x3e, the coordinates which locate the masses are obtained from Eq. (3), i.e. r1e = x1e and r2e = x1e + x3e + L. Next, a static nonlinear feedback control laws to drive the masses to the equilibrium point is going to be analyzed. Let us suppose that the values of A1 and xd in Eq. (19) are chosen to obtain chaotic behaviour. Since the attractor is a dense set, there will be a trajectory intersecting a rectangle X defined as follows: X ¼ fjx3 ðtÞ  x3e j 6 ra1 ; jx4 ðtÞj 6 ra2 g

ð20Þ

where ra1 and ra2 are previously chosen small values. At this moment, it will be possible to destroy the chaotic motion by commuting to a new control law defined by: ) u1 ðtÞ ¼ as K q x33 ðtÞ  B1 x4 ðtÞ þ F ð21Þ e u2 ðtÞ ¼ K p x3 ðtÞ  B2 x4 ðtÞ þ 2m F m1 where F may be considered to be a constant force applied by the actuators of the device (see Fig. 1a). Note that, due to the particular coupling of the oscillators, the terms as K q x33 ðtÞ  B1 x4 ðtÞ and Kpx3(t)  B2x4(t) are necessary to assure the compatibility of Eq. (7) at the equilibrium point. So, the control law (21) should be considered as simple as possible. To analyze the system defined by Eq. (7) with the control law (21), it is convenient to research the behaviour of the oscillators associated to the state variables x1(t), x2(t) and x3(t), x4(t) separately. Substituting Eq. (21) into Eq. (7), the oscillator associated to the state variables x3(t) and x4(t) can be written as follows: ) x_ 3 ðtÞ ¼ x4 ðtÞ Kp K K Kq B B1 B2 a2 ¼  ; a3 ¼  ; b2 ¼   ð22Þ x_ 4 ðtÞ ¼ a2 x3 ðtÞ  as a3 x33 ðtÞ  b2 x4 ðtÞ þ mF1 m1 me me m1 me m1 me From Eq. (22), the characteristic equation which defines the jacobian matrix at the equilibrium point x3e gives the following possibilities of reachable and unreachable equilibrium points: ( saddle points if a2 > 3as a3 x23e 2 2 ð23Þ k þ b2 k  ða2  3as a3 x3e Þ ¼ 0 ) focus or stable node if a2 < 3as a3 x23e Consequently, there is an interval of unreachable equilibrium points around the origin, which is given by: rffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffi a2 a2  < x3e < þ 3as a3 3as a3

ð24Þ

On the other hand, from Eq. (22) it is possible to represent the variation of the control force F as a function of the equilibrium point x3e, i.e: F ¼ m1 as a3 x33e  m1 a2 x3e

ð25Þ

Taking into account the parameter values indicated in legend of Fig. 1, the plot of Eq. (25) is shown in Fig. 3. Note that there is a zone of unreachable points (±17.563 · 106), a zone with three equilibrium points and two symmetrical zones around the origin with one equilibrium point. Therefore, choosing a set point x3e, the value of F can be calculated from Eq. (25). However, if the value of F is inside the zone Fmin < F < Fmax, a new equilibrium point with different sign may appear. In this case, if the state variables are in the basin attraction of this new point, just when the control law is commuted the system will be driven to this point, and the control law (21) will not drive the masses to the desired set point. This difficulty may be overcome by choosing equilibrium points outside of this zone. To analyze the control law of Eq. (21), the system has been simulated twice, as shown at Fig. 4. First, the system is driven to chaotic behaviour by means of the control law (19). Next, at t = 0.5 s, the control law (19) is commuted to the control law (21). Fig. 4a shows that the set point defined by the values x1e = 45.678 · 106 m and x3e = 18.999 · 106 m

1362

Manuel F. Pe´rez Polo et al. / Chaos, Solitons and Fractals 39 (2009) 1356–1370

Fig. 3. Plot of the reachable and unreachable equilibrium points from Eq. (25), i.e. values of the constant force F vs. the equilibrium point x3e.

Fig. 4. Simulation results from the control law (21). Parameter values are indicated in legend of Fig. 1. Set point: r1e = x1e = 45.678 · 106 m, x3e = 18.999 · 106 m, r2e = x3e + r1e + L = 74.677 · 106 m. (a) The commutation of the control law (19)–(21) is forced to be at t = 0.5 s with d = 500. The equilibrium points x1er and x3er are not coincident with the set points x1e and x3e. (b) Strange attractor of case (a). (c) The commutation of the control law (21) is forced to be at t = 0.6 s with d = 500. The equilibrium point is reached. (d) Strange attractor of case (c).

is not reached. This is due to the fact that the value of the force is F = 0.0242N, which is. inside the zone Fmin < F < Fmax (see Fig. 3). Consequently, another stable equilibrium point appears at x3er = 34.362 · 106 m. Fig. 4b shows the strange attractor and how the desired set point is not reached because of the presence of x3er. Fig. 4c and d show the same case, but now the control is commuted at t = 0.6 s, when the system is in the attraction basin of the set point. Fig. 5 shows two cases in which the control law (19) is commuted to the control law (21) when a chaotic trajectory intersects the capture zone X defined in Eq. (20). In Fig. 5a the capture zone around the equilibrium point x3e is very

Manuel F. Pe´rez Polo et al. / Chaos, Solitons and Fractals 39 (2009) 1356–1370

1363

Fig. 5. Simulation results from the control law (21). Parameter values are indicated in legend of Figs. 1 and 4. The commutation of the control law (19)–(21) is fulfilled when intersecting the capture zone (Eq. (20)) with d = 500. (a) Capture zone: ra1 = 3 · 107 m, ra2 = 3 · 105 m/s. The equilibrium point is reached. (b) Control forces of case (a). (c) Capture zone: ra1 = 5 · 106 m, ra2 = 3 · 105 m/s. The equilibrium point is not reached (d) Strange attractor of case (c).

small and therefore the set point is reached. In Fig. 5b, the control forces defined by Eq. (8) are plotted, showing that they are in the magnitude order of the data reported in the bibliography. Since the size of the capture zone is not adequate, the same problem of Fig. 4a and b appears in Fig. 5c and d, i.e. the desired set point is not reached.

5. Positioning of masses with harmonic base excitation The purpose of this section is to change the control laws of the previous section when harmonic platform excitations appear in the device. The start point must be Eq. (7), where the external harmonic disturbances are defined by the following expressions: ) x1 ðtÞh_ 2 ¼ x1 ðtÞA2h cos2 xh t ð26Þ ½x3 ðtÞ þ Lh_ 2 ¼ ½x3 ðtÞ þ LA2 cos2 xh t h

Let us suppose that the amplitude Ah and frequency xh can be measured by means of adequate sensors, such as accelerometers type MEMS integrated in the same device. Then, it may be possible to construct the harmonic disturbance as a function of time and inject it into the device to cancel the undesired platform vibration. This compensation procedure may be considered as a feedforward technique [25], which works properly if the disturbance can be measured with precision. In this work we assume that this is the case. The first problem consists of determining how to change the chaotic control law defined by Eq. (19) to obtain chaotic behaviour, even with external harmonic disturbances. To do this end, note that the mean value of signal h_ 2 ¼ A2h cos2 xh t is not zero, i.e: Z T A2 x2 _2 1 h ¼ A2h cos2 xh t dt ¼ h h ð27Þ 2 T 0 Consequently, the terms x1 ðtÞA2h x2h =2; ½x3 ðtÞ þ LA2h x2h =2 appear in the second and four Eq. (7) respectively. This mean value gives a continuous shift of masses m1 and m2, making the system unstable. In this situation, the control law (19) cannot drive the masses to chaotic behaviour. To avoid this problem, the control signals are modified to compensate the mean value (27), i.e.:

Manuel F. Pe´rez Polo et al. / Chaos, Solitons and Fractals 39 (2009) 1356–1370

1364

u1 ðtÞ ¼ as K q x33 ðtÞ  B1 x4 ðtÞ  m1 x1 ðtÞ u2 ðtÞ ¼ K p x3 ðtÞ  B2 x4 ðtÞ  me

A2h x2h 2

A2h x2h 2

9 =

þ A1 cos xd t

e ½x1 ðtÞ þ x3 ðtÞ þ L þ 2m A1 cos xd t ; m1

ð28Þ

By means of Eq. (28), the DC (Eq. (27)) components of x1 ðtÞh_ 2 ; ½x3 ðtÞ þ Lh_ 2 are removed, and if the harmonic base disturbance exists, it will be added to the harmonic vibration which produces chaotic behaviour, without shift of the masses. The same reasoning can be used respect to the control law (19), i.e. 9 A2 x2 = u1 ðtÞ ¼ as K q x33 ðtÞ  B1 x4 ðtÞ  m1 x1 ðtÞ h2 h þ F ð29Þ 2 2 A x u ðtÞ ¼ K x ðtÞ  B x ðtÞ  m h h ½x ðtÞ þ x ðtÞ þ L þ 2me F ; 2

p 3

2 4

e

2

1

3

m1

From control laws (28) and (29) only the DC component of harmonic base disturbances is compensated. So to remove the harmonic vibration completely, new harmonic components are necessary. In this case, Eqs. (19) and (21) must be modifies by adding the terms: ) fd1 ðtÞ ¼ m1 x1 ðtÞA2h x2h sin2 xh t ð30Þ fd2 ðtÞ ¼ me ½x1 ðtÞ þ x3 ðtÞ þ LA2h x2h sin2 xh t Figs. 6 and 7 show the simulation results when there is only one equilibrium point, i.e. F > Fmin; F > Fmax (see Fig. 3). Fig. 6a shows chaotic oscillations of the masses since the control law (21) is being applied without harmonic base excitation. At t  0.75 s, a chaotic trajectory intersects the capture zone, and the control law (30) drives the system to the set point. At t = 0.9 s, a harmonic base disturbance appears, and the state variables x1(t) and x3(t) remain oscillating. The presence of a small limit cycle in the phase plane x3(t)  x4(t) is shown in Fig. 6b, which disappears when the control is applied and eventually the set point is reached. At t = 1.1 s, the control laws (28)–(30) are applied, compensating the harmonic disturbances and driving the system to the set point, as shown in Fig. 7a. Note the appearance of a sensitive dependence zone, which disappears when the chaotic motion is destroyed as shown in Fig. 7b. Fig. 7c shows

Fig. 6. Simulation results from the control law (29) and (30). Parameter values are indicated in legend of Fig. 1. Set point: r1e = x1e = 25.678 · 106 m, x3e = 48.999 · 106 m, r2e = x3e + r1e + L = 84.677 · 106 m. The commutation of the control law (28) to (29) and (30) is fulfilled when intersecting the capture zone (Eq. (20)) with d = 500. (a) Capture zone: ra1 = 2 · 107 m, ra2 = 1 · 103 m/s. Harmonic disturbance with Ah = 1.1 rad and xh = 500 rad/s is applied at t = 0.9 s and compensated at t = 1.1 s. The equilibrium point is reached (b) Strange attractor of case (a) and limit cycle due to the disturbance.

Manuel F. Pe´rez Polo et al. / Chaos, Solitons and Fractals 39 (2009) 1356–1370

1365

Fig. 7. Simulation results from Fig. 6. (a) Chaotic oscillations of variable x1(t), steady state, harmonic disturbances and cancellation of the harmonic disturbances. (b) Sensitive dependence of case (a). (c) Control signal u1(t) of mass m1 and oscillating control due to the harmonic base disturbances. (d) Variation of the Lyapunov exponents vs. time.

that the value of the control signal u1(t) is between ±0.5 N, which can be considered to be adequate in these MEMS devices. Fig. 7d shows that there is a positive Lyapunov exponent, but all Lyapunov exponents are eventually negative when the chaotic motion is suppressed.

6. Geometric nonlinear control The control laws deduced in Sections 4 and 5 are not systematic, i.e. they are obtained from a physical point of view. However, this section focuses on developing a systematic approach to the positioning problem by means of nonlinear control techniques. The starting point will be the system equations written in affine form [20,21] i.e: 9 m P x_ ðtÞ ¼ f ½xðtÞ þ gi ½xðtÞui ðtÞ = ð31Þ i¼1 ; y½xðtÞ ¼ h½xðtÞ where x(t) is the state vector, f[x(t)] and gi[x(t)] are vector fields and y[x(t)] are the output nonlinear signals defined by the designer. Then, Eq. (7) in affine form can be written as follows: 9 3 2 3 2 3 3 2 2 x2 ðtÞ 0 0 > x_ 1 ðtÞ > > > 7 6 1 7 K K 3 B > 7 6 x_ ðtÞ 7 6 6 x ðtÞ þ a x ðtÞ þ x ðtÞ > 0 7 6 3 s 4 > m1 3 m1 7 6 2 7 6 m1 6 m1 7 6 > 7 þ u u ðtÞ þ ðtÞ > 7¼6 7 7 6 6 6 1 1 = 7 0 5 5 4 x_ 3 ðtÞ 5 4 4 4 0 x4 ðtÞ 5 ð32Þ 1 1 K K 3 B >  m1  me  me x3 ðtÞ  as me x3 ðtÞ þ me x4 ðtÞ x_ 4 ðtÞ > > >



> > > h1 ½xðtÞ x1 ðtÞ > > y½xðtÞ ¼ h½xðtÞ ¼ ¼ ; h2 ½xðtÞ x3 ðtÞ where it is assumed that there is no base excitation. Note that the output signals have been defined from the state variables x1(t) and x3(t). From Eqs. (31) and (32), it is easy to identify the vector fields f[x(t)] and gi[x(t)]. The control laws

Manuel F. Pe´rez Polo et al. / Chaos, Solitons and Fractals 39 (2009) 1356–1370

1366

are obtained from the exact input–output linearization. For this purpose, taking the derivative respect to time in the second Eq. (32), it is deduced that: _ ¼ yðtÞ

oy dxðtÞ ohðxÞ ¼ ff ½xðtÞ þ g1 ½xðtÞu1 ðtÞ þ g2 ½xðtÞu2 ðtÞg ¼ Lf hðxÞ þ Lg1 hðxÞu1 ðtÞ þ Lg2 hðxÞu2 ðtÞ ox dt ox

ð33Þ

where the Lie derivative nomenclature has been used, i.e. Lfh(x) = (oh(x)/ox)f(x) [20,21]. From Eqs. (31) and (32), the Lie derivatives can be written as follows: 3 2 0 6 1=m 7 oh1 ðxÞ 1 7 6 ð34Þ Lg1 h1 ðxÞ ¼ g1 ðxÞ ¼ ½ 1 0 0 0 6 7¼0 4 0 5 ox 1=m1 Similarly, by direct computation it is obtained that: Lg2 h1 ðxÞ ¼ ½ 1 0 0 0 g2 ðxÞ ¼ 0 Lg1 h2 ðxÞ ¼ ½ 0 0 1 0 g1 ðxÞ ¼ 0;

Lg1 h2 ðxÞ ¼ ½ 0 0 1 0 g1 ðxÞ ¼ 0

ð35Þ

Since all Lie derivatives are zero, it is necessary to take the derivative respect to time once again, i.e.: €y 1 ðtÞ ¼ L2f hðxÞ þ Lg1 Lf h1 ðxÞu1 ðtÞ þ Lg2 h1 ðxÞu2 ðtÞ

ð36Þ

oLf h1 ðxÞ Lg1 Lf h1 ðxÞ ¼ g1 ðxÞ ¼ 1=m1 ; ox

ð37Þ

oLf h1 ðxÞ Lg2 Lf h1 ðxÞ ¼ g2 ðxÞ ¼ 0 ox

Taking into account that the same calculations can be realized from the output y2(t), the following matrix equation is obtained: #



" 2 Lf h1 ðxÞ €y 1 ðtÞ 1=m1 0 u1 ðtÞ þ ¼ ð38Þ €y 2 ðtÞ 1=m1 1=me u2 ðtÞ L2f h2 ðxÞ where the values of the Lie derivatives are the following: ) L2f h1 ðxÞ ¼ mK1 x3 ðtÞ þ as mK1 x33 ðtÞ þ mB1 x4 ðtÞ L2f h2 ðxÞ ¼  mKe x3 ðtÞ  as mKe x33 ðtÞ  mBe x4 ðtÞ From Eqs. (38) and (39), the input signals as a function of the output derivatives can be written as follows: 9

= u1 ðtÞ ¼ m1 €y 1 ðtÞ  Kx3 ðtÞ þ as Kx33 ðtÞ þ Bx4 ðtÞ  

u2 ðtÞ ¼ me ½€y 1 ðtÞ þ €y 2 ðtÞ þ 1  mm1e Kx3 ðtÞ þ as Kx33 ðtÞ þ Bx4 ðtÞ ;

ð39Þ

ð40Þ

On the other hand, introducing the auxiliary variables: z1 ðtÞ ¼ y 1 ðtÞ;

z2 ðtÞ ¼ y_ 1 ðtÞ;

z3 ðtÞ ¼ y 2 ðtÞ;

the following linear system is deduced: 3 2 2 0 z_ 1 ðtÞ 6 z_ ðtÞ 7 6 0 6 2 7 6 z_ ðtÞ ¼ AzðtÞ þ B€y ðtÞ ) 6 7¼6 4 z_ 3 ðtÞ 5 4 0 0 z_ 4 ðtÞ

1 0 0 0

z4 ðtÞ ¼ y_ 2 ðtÞ

0 0 0 0

32 3 3 2 z1 ðtÞ 0 0 0 7

7 6 6 07 76 z2 ðtÞ 7 6 1 0 7 €y 1 ðtÞ 76 7 7þ6 1 54 z3 ðtÞ 5 4 0 0 5 €y 2 ðtÞ 0 0 1 z4 ðtÞ

The output equations of system (42) are chosen as follows: 3 2 z1 ðtÞ



6 1 0 0 0 6 z2 ðtÞ 7 y 1 ðtÞ 7 ¼ yðtÞ ¼ CzðtÞ ) 7 6 y 2 ðtÞ 0 0 1 0 4 z3 ðtÞ 5

ð41Þ

ð42Þ

ð43Þ

z4 ðtÞ Since the output derivatives y_ 1 ðtÞ, €y 2 ðtÞ are independents, it is possible to apply the pole-placement technique [25,26] taking into account a linear feedback law as follows:

Manuel F. Pe´rez Polo et al. / Chaos, Solitons and Fractals 39 (2009) 1356–1370

€y ðtÞ ¼ KzðtÞ )





€y 1 ðtÞ K 11 ¼ €y 2 ðtÞ K 21

K 12 K 22

K 13 K 23

3 2 z1 ðtÞ 6 K 14 6 z2 ðtÞ 7 7 7 6 K 24 4 z3 ðtÞ 5

1367

ð44Þ

z4 ðtÞ Substituting Eq. (44) into Eq. (42), it is obtained a system which the constant coefficientsKij are chosen to obtain a Hurwitz type matrix. From the theory it is known that the matrix K in Eq. (44) can always be determined if

of linear systems .. .. 2 .. 3 the matrix M defined as: M ¼ B.AB.A B.A B has rank n = 4. The constants K can be determined by applying the ij

Ackermann’s formula for multivariable systems. Details can be found in references [25,26]. From the previous considerations, the derivatives of the output signals can be written as follows:  €y 1 ðtÞ ¼ K 11 ½x1 ðtÞ  x1e   K 12 x2 ðtÞ  K 13 ½x3 ðtÞ  x3e   K 14 x4 ðtÞ ð45Þ €y 2 ðtÞ ¼ K 21 ½x1 ðtÞ  x1e   K 22 x2 ðtÞ  K 23 ½x3 ðtÞ  x3e   K 24 x4 ðtÞ where deviation variables have been introduced. Consequently, from Eqs. (40) and (45), the control law is obtained. It is interesting to remark that if rank(M) = 4 and the derivatives of the outputs are equal to the number of state variables, the previous method always gives a systematic procedure to obtain the control signals. In the absence of external harmonic disturbances, by applying Eq. (19), it may be possible to obtain chaotic behaviour. Then, when a chaotic trajectory intersects a capture zone, by commutating to the control laws (40) and (45), it will be possible to drive the masses to the equilibrium point. The appearance of harmonic base excitations can be considered as in the previous section, i.e. taking into account the DC value of the disturbances (Eq. (27)) and the possibility of removing them by means of cancellation. For example, to remove the the harmonic base disturbance, Eq. (40) must been changed by the new ones: 9 = u1 ðtÞ ¼ m1 €y 1 ðtÞ  ½Kx3 ðtÞ þ as Kx33 ðtÞ þ Bx4 ðtÞ  m1 x1 ðtÞA2h x2h sin2 xh t   ð46Þ u2 ðtÞ ¼ me ½€y 1 ðtÞ þ €y 2 ðtÞ þ 1  mm1e ½Kx3 ðtÞ þ as Kx33 ðtÞ þ Bx4 ðtÞ  me A2h x2h ½x1 ðtÞ þ x3 ðtÞ þ L sin2 xh t ; The dynamical behaviour of the masses defined by Eq. (7) with the previous control laws is shown in Fig. 8. Fig. 8a shows that the control law (19) gives chaotic behaviour. At approximately 0.54 s, a chaotic trajectory enters in the capture zone and the commutation to the control laws (41) and (46) take places driving the system to the set point even with harmonic base excitation. In Fig. 8b and d, the set point is outside the strange attractors, while in Fig. 8c, the set point is inside them. The eigenvalues assigned to calculate the constants K of Eq. (44) are ½ 100  300i 500 600 , and the corresponding values are:

2:7973 0:0071 7:5813 0:0063 K ¼ 105 9:4079 0:0094 8:3380 0:0189 The simulation results show that the values of the control forces and accelerations are similar to the previous one.

7. Comparison of the control laws at very short times The simulation results presented in the previous sections show that the researched control laws can drive the MEMS device with nanometre precision. The simulation time was approximately 1 s, in order to visualize how the control laws work in different situations, i.e. depending on the presence or absence of harmonic base disturbances. In practical cases, however, these devices will be designed to operate at very short times, i.e. the read–write time of the probe tips will probably be of milliseconds. To research this issue, in this section we are going to compare some of the previously deduced control laws at very short times. The idea is to avoid an excessive seek time of the read–write heads as much as possible since, depending on the applications this may be a disadvantage. This problem can be overcome by increasing the capture zone to decrease the seek time. We are going to consider the more relevant issues when the control laws (19) or (28) are commuted to the control laws (21) and (30), depending on whether the harmonic base excitation is considered or not. Fig. 9 shows the results obtained with the control law (21) without taking into account harmonic base excitations. As we saw in Section 4, the disadvantage of this control lies in the existence of an unreachable zone in the magnetic media. The set point is quickly reached by changing the parameter d defined in Eq. (17) to values between 1000 and 2000. Since the initial conditions are appropriate, mass m1 reaches the set point with a positioning time of 7.5 ms, as shown in Fig. 9a. The commutation of the control law is clearly shown in Fig. 9b. Fig. 9c and d show the change of the forces

1368

Manuel F. Pe´rez Polo et al. / Chaos, Solitons and Fractals 39 (2009) 1356–1370

Fig. 8. Simulation results from the control law (48) and (45). Parameters values are indicated in legend of Fig. 1. Set point: r1e = x1e = 77.899 · 106 m, x3e = 15.432 · 106 m, r2e = x3e + r1e + L = 72.467 · 106 m. The commutation of the control law (19) to (40) and (45) are fulfilled when a chaotic trajectory intersecting the capture zone (Eq. (20)) with d = 500. (a) Capture zone: ra1 = 5 · 107 m, ra2 = 2 · 103 m/s. Harmonic disturbance with Ah = 1.3 rad, xh = 600 rad/s is applied at t = 0.8 s and compensated at t = 0.9 s. The equilibrium point is reached (b)–(d) Strange attractors of case (a) and limit cycles due to the disturbance.

Fig. 9. Simulation results from the control law (21). Parameters values are indicated in legend of Figs. Figs. 1 and 4. Set point: r1e = x1e = 40.987 · 106 m, x3e = 32.345 · 106 m, r2e = x3e + r1e + L = 83.332 · 106 m. Simulation time = 0.01 s, simulation step = 1 · 106 s. The commutation of the control law (19)–(21) is fulfilled when a chaotic trajectory intersects the capture zone (Eq. (20)) with d = 2000. Capture zone. ra1 = 1 · 105 m, ra2 = 0.05 m/s. There are not harmonic disturbances. (a) Simulation values of x1(t) and x3(t). (b) Abrupt commutation of the control law. (c) Control forces F1(t) and F2(t). (d) Accelerations of masses m1 and m2.

Manuel F. Pe´rez Polo et al. / Chaos, Solitons and Fractals 39 (2009) 1356–1370

1369

and accelerations when the commutation from control law (19)–(29) has been fulfilled. Note that the length of the chaotic trajectory is very small because the capture zone is very large. Similar results can be obtained from the control laws (40) and (45). Therefore by choosing and adequate capture zone, it will be possible positioning the masses with settling times of milliseconds and nanometre precision.

8. Conclusions This paper develops a nonlinear mechanical model for MEMS-based storage devices in an attempt to obtain chaotic behaviour and nonlinear control in a unified framework. The model is formed by a damping and a hard spring which interconnects two masses. The proposal is to design an adequate control law to obtain positioning with nanometre precision and settling times of milliseconds, as well as typical velocities and accelerations of actual MEMS devices. By using the adequate state variables, it is shown that the model can be represented by two asymmetrically coupled oscillators. The paper shows that, as a result of this particular connection, two Duffing’s oscillators can be obtained. Next, by adding external periodic forces, and from the Melnikov method, the sufficient condition for chaotic motion is applied to determine a range for the amplitude and frequency of the external forces, which give chaotic motion. The chaotic motion of the oscillators is verified by means of the appearance of strange attractors, power spectral density and the computation of Lyapunov exponents. The simulation results are indirectly verified by means of the sum of the Lyapunov exponents, which approximately equals the divergence of the vector system field. The paper shows that the size of the strange attractor, the velocities, accelerations and the control forces are in accordance with the values of actual MEMS. The chaotic motion is used to reduce the velocities, accelerations and forces by commuting to an adequate control law when the masses are very close to the set point. To do this end, a small capture zone is defined. This zone will always be intersected by a chaotic trajectory, and then, the commutation law will arise. So, it will be possible to solve the positioning problem with an adequate control effort and very short seek times. The commutation from chaotic to regular motion has been researched considering two control laws: (i) A static nonlinear control law with constant forces and (ii) Exact input–output linearization by means of nonlinear control based on differential geometry. Case (i) gives the most simple control law, but with the drawbacks of the zones, which are unreachable for the read–write heads. The control law in case (i) is little systematic, so a new geometric control law in connection with the pole-placement technique has been deduced. The compensation of the harmonic base excitation is realized by using the same control laws. Assuming that the amplitude and frequency of the harmonic disturbances can be measured, a compensating feedforward law can be added to the previous control laws to cancel the harmonic disturbances. The simulation results show that, by means of the previously designed control laws, it is possible to cancel the harmonic disturbances and obtain the desired set point, with an appropriate dynamic behaviour. Finally, the paper shows that the previous control laws can also be applied at very short times if the capture zone is large. Then, when a short time has elapsed, a chaotic trajectory will intersect the capture zone and the commutation of the control law will destroy the chaotic motion to drive the masses to the desired set point. The simulation results show that a good dynamical behaviour, with minimum control effort and settling times of ms is obtained. This paper can be extended in many directions. For example, the proposed mathematical model can be modified by adding the effect of nonlinear anchor springs, or by considering the elastic nonlinear hereditary behaviour of the mass sled. In this case, new nonlinear models arise, and there exists a possibility of researching new chaotic behaviours in connection with new control laws, which could be designed from the ideas presented in this paper.

References [1] [2] [3] [4] [5]

El Naschie MS. Chaos and fractals in nano and quantum technology. Chaos Solitons & Fractals 1998;9(10):1793–802. El Naschie MS. Nanotechnology for the developing world. Chaos Solitons & Fractals 2006;30:769–73. Rosen J. Machining in the micro domain. Mech Eng 1989;111(3):40–6. O’Connor L. MEMS microelectromechanical systems. Mech Eng 1992;114(2):40–7. Reiley TC, Fan LS, Mamin HJ. Micromechanical structures for data storage. Journal of Microelectronic Engineering 1995;27:495–8. [6] Griffin JL, Schlosser SW, Ganger GR, Nagle DF. Modelling and performance of MEMS-based storage devices. In: Proceedings of ACM Sigmetrics 2000, Santa Clara, California, June, 2000. [7] Maithirpala DMS, Kawade BD, Berg JM. General modeling and control framework for electrostatically actuated mechanical systems. Int J Robust Nonlinear Control 2005:1–19.

1370

Manuel F. Pe´rez Polo et al. / Chaos, Solitons and Fractals 39 (2009) 1356–1370

[8] Liu S, Davidson A, Lin Q. Simulation studies on nonlinear dynamics and chaos in a MEMS cantilever. J Micromech Microeng 2004;14:1064–73. [9] Le´vesque L. Revisiting the coupled-mass system and analogy with a simple band gap structure. Eur J Phys 2006:133–45. [10] Sinha SC, Gourdon E, Zhang Y. Control of time-periodic systems via symbolic computation with application to chaos control. Commun Nonlinear Sci Numer Simulat 2005;10:835–54. [11] Haitao M, Deshmukh V, Butcher E, Averina V. Delayed state feedback and chaos control for time-periodic systems via symbolic approach. Commun Nonlinear Sci Numer Simulat 2005;10:479–97. [12] Wiggins S. Introduction to applied nonlinear dynamical systems and chaos. New York: Springer-Verlag; 2003. p. 687–720. [13] Guckenheimer J, Holmes P. Nonlinear oscillations, dynamical systems and bifurcations of vector fields. New York: SpringerVerlag; 1983. p. 173–7, 191–3, 198–201. [14] Lichtenberg AJ, Lieberman MA. Regular and chaotic dynamics. New York: Springer Verlag; 1992. p. 296–312, 315–20, 560–9. [15] Hilborn RC. Chaos and nonlinear dynamics. Oxford: Oxford University Press; 2000. p. 579–83. [16] Marsden JE, Ratiu TS. Introduction to mechanics and symmetry. New York: Springer Verlag; 1999. p. 94–103. [17] Lu Michal SC, Fedder GK. Position control of parallel-plate microactuators for probe-based data storage. IEEE J Microelectromech Syst 2004;13(5):759–69. [18] Fedder GK. Multimode digital control of a suspended polysilicon microstructure. IEEE J Microelectromech Syst 1996;5(4):283–97. [19] de Sudipto K, Aluru NR. Complex nonlinear oscillations in electrostatically actuated microstructures. IEEE J Microelectromech Syst 2006;5(2):355–69. [20] Nejmeijer H, Van der Schaft. Nonlinear control. London: Springer-Verlag; 1990. [21] Isidori A. Nonlinear control system. London: Springer-Verlag; 1995. p. 137–89. [22] Goldstein H. Classical mechanics. New York: Addison Wesley; 1980. [23] Benettin G, Galgani L, Giorgilly A, Strelcyn JM. Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them, Part I: Theory. Meccanica 1980;15:9–20. [24] Benettin G, Galgani L, Giorgilly A, Strelcyn JM. Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them, Part II: numerical applications. Meccanica 1980;15:21–30. [25] Albertos P, Sala A. Multivariable control systems. London: Springer; 2004. p. 165–88. [26] Antsaklis PJ, Michel AN. Linear systems. New York: McGraw-Hill; 1997. p. 226–47, 326–40.