Chaos, Solitons and Fractals 12 (2001) 49±66
www.elsevier.nl/locate/chaos
Controlling the nonlinear dynamics of a beam system M.F. Heertjes *, M.J.G. Van De Molengraft Faculty of Mechanical Engineering, Eindhoven University of Technology, W-Hoog 1.143, P.O. Box 513, Den Dolech 2, 5600 MB Eindhoven, Netherlands Accepted 23 August 1999
Abstract Control based on linear error feedback is applied to reduce vibration amplitudes in a piecewise linear beam system. Hereto small amplitude 1-periodic solutions are stabilized wherever they coexist with two or more long-term solutions. In theory, no control eort is required to maintain the 1-periodic response once it has been stabilized. For the beam system, 1-periodic solutions are stabilized by feedback at one location along the beam. Feedback is represented by servo-stiness or servo-damping which results from increasing two corresponding control parameters. At appropriate levels of these parameters local, or global, asymptotic stability (of the zeroequilibrium) of the error dynamics, i.e. stability of the underlying 1-periodic solutions, can be guaranteed. Local asymptotic stability can be guaranteed for a large range of actuator locations and excitation frequencies and is indicated by bifurcations. Global asymptotic stability can only be guaranteed for a limited range of actuator locations on the basis of the well-known circle criterion. The dierence between local and global asymptotic stability in terms of the required values for the control parameters can be signi®cant, and may result in large dierences in control performance. Ó 2000 Elsevier Science Ltd. All rights reserved.
1. Introduction Periodically excited nonlinear systems, such as suspension bridges [5,14], gear-boxes [12,17], or atomic force microscopes [1,22], appear frequently in engineering practice. Such systems are known to possess two or more coexisting long-term solutions at certain regions of the system parameter space. One of these solutions usually has preferable features, for example, smaller vibration amplitudes, that can prevent the system from damage or wear. This favors the idea of choosing one of the coexisting solutions as the desired trajectory for control as is common in the control of chaos [19]. Especially, because the choice of a natural solution as the desired trajectory for control may result in limited control eort. Hence no further control force is needed once the desired long-term response has been stabilized. Stabilizing coexisting long-term solutions has received much attention in the research on chaos [3]. Mostly, this is due to the ergodic property which assures that almost all trajectories of the system will reach in ®nite time any small neighborhood of any point in the chaotic attractor. However, the control of chaos suers from two major de®ciencies with respect to vibration amplitude reduction. Firstly, it is generally undesirable to exploit the ergodic property when preventing damage or wear because, in this case, the exposure to harmful vibrations must be minimized. Secondly, harmful long-term vibrations are often not chaotic. To overcome these de®ciencies, a control design is required that may be applied globally rather than the commonly used designs based upon linearization that can only be applied locally. To this extent, a control design based on linear error feedback will be applied to a piecewise linear beam system.
*
Corresponding author. Fax: +31-40-2461418. E-mail address:
[email protected] (M.F. Heertjes).
0960-0779/00/$ - see front matter Ó 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 9 6 0 - 0 7 7 9 ( 9 9 ) 0 0 1 6 4 - 2
50
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66
The piecewise linear beam system under study consists of a supported beam which is excited harmonically, and which is allowed to impact at its middle on a one-sided spring [18]. This system is related to the so-called impact oscillator [2] but diers at some major points. Firstly, the beam will be modeled by a 3degree-of-freedom (DOF) model whereas the impact oscillator mostly refers to a single-DOF model. Secondly, the one-sided spring has ®nite stiness whereas the impact oscillator impacts on an in®nitely sti wall. The piecewise linear beam system shows coexisting long-term solutions. Especially, 1-periodic solutions are of interest as they exist throughout the considered frequency-range and have the smallest vibration amplitude compared to the other coexisting long-term solutions [10]. This, together with the expected convergence of the control force to zero, makes 1-periodic solutions desirable trajectories to track. Linear error feedback, which will be used to stabilize 1-periodic solutions, will be applied at only one location along the beam. At this location, the error dynamics with respect to 1-periodic solutions can be changed by increasing two control parameters: one for the amount of servo-damping, and one for the amount of servo-stiness. Under certain conditions, the error dynamics are locally, or globally, asymptotically stable, i.e. the underlying 1-periodic solutions are locally, or globally, asymptotically stable. Global asymptotic stability may be guaranteed on the basis of the circle criterion [13,15]. Unfortunately, the circle criterion only provides a sucient condition for stability which for the beam system applies to a limited range of actuator locations. Moreover, to guarantee global asymptotic stability, large values for the control parameters may be required. For much smaller values of the control parameters bifurcations occur that indicate the disappearance of coexisting solutions. These values can be found by local stability analysis and may provide signi®cantly better results in terms of control eort and settling time of the error dynamics. Besides, local asymptotic stability can be guaranteed for a larger range of actuator locations and excitation frequencies. 2. Piecewise linear beam system The lab-scale setup of the piecewise linear beam system is depicted in Fig. 1 and contains four characteristic entities: a beam, a one-sided spring, an excitation force v, and a control force u. The beam is an elastic structure supported by leaf springs at both ends. The one-sided spring is constructed as a second beam that adds stiness to the system for positive displacements of the middle of the main beam [11,18]. The excitation force v is a sine-shaped force acting on the middle of the beam and is generated by a massunbalance. The mass-unbalance is driven by a tacho-controlled motor at a user-de®ned rotation speed. The nonlinear behavior of the beam system can be in¯uenced by a control force u that is applied somewhere beam length = 1.3 [m] 0.65 [m] L=0.382[m] W=0.006[m] H=0.02 [m]
0.325 [m]
a
one-sided spring
leaf spring L=0.09[m] W=0.001[m] H=0.075 [m] 2
nodes 3-33
34
nodes 35-66
L=1.3[m] W=0.01[m] H=0.1 [m] 67
beam driving shaft
qact
qmid
u
v actuator dSPACE/ simulink
y
mass-unbalance a
motor cross-section: a
z x
LVDT
beam
y mass-unbalance force transducer
z
x
accelerometer
one-sided spring
Fig. 1. Graphical representation of the piecewise linear beam system.
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66
51
along the beam. The control force is generated by an actuator represented by a shaker±ampli®er combination [9]. Both the control force u and the excitation force v are measured by piezo-electric force transducers with an accuracy of 10ÿ1 N. Furthermore, the displacements qact at the actuator location and qmid at the middle of the beam are measured by linear variable dierential transformers (LVDTs), whereas the qmid are measured with piezo-electric accelerometers. The accuracy of corresponding accelerations qact and these measurements are 10ÿ6 m and 10ÿ2 msÿ2 , respectively. Acquisition and processing of the measured data together with the on-line computation of the control force are performed with the real-time hardware/software environment dSPACE [6] combined with Matlab's Simulink [20]. For control design, a simple though accurate model of the beam system is required. An accurate model of the beam is built that contains 152-DOF. This model is reduced by a component mode synthesis technique [4] to reduce the required computational eort in future control experiments. The loss of accuracy in the reduced model mainly applies to the high frequency behavior of the linear beam model. It is assumed that this behavior has little in¯uence on the low frequency response of the nonlinear beam system. A 3-DOF model of the beam satis®es this assumption in the frequency-range of interest. The 3-DOF model contains two physical DOFs and one non-physical DOF. The physical DOFs, qact and qmid , are needed to incorporate the external forces into the nonlinear model of the beam system. The external forces are: the excitation force, the one-sided spring force, and the control force. The non-physical DOF n is needed for an accurate description of the ®rst eigenmode of the linear beam [7]. The nonlinear model of the beam system is given by the following equations of motion Mq Bq_ Kq knl
qmid h2 hT2 q h2 v h1 u; with
qmid
1 0
if qmid > 0; if qmid 6 0;
1
2
and 2
v v
t c
2pfexc cos
2pfexc t; T
3
where q qact qmid qn , M, B, and K are positive de®nite 3 3 matrices with B based on modal damping f 0:05, which is chosen equal for all DOFs, hi is a distribution matrix with zeros except for the ith entry which equals one, knl represents the stiness of the one-sided spring, c an excitation constant, fexc the excitation frequency, and u the input; see Appendix A for numerical values of a 3-DOF model. The validity of this model is assessed in the time-domain by comparing numerical with experimental results. A time-domain validation at a ®xed excitation frequency of 35 Hz is shown in Fig. 2, where time-series for a long-term 2-periodic solution at an excitation frequency of 35 Hz are depicted; the actuator was located at a quarter of the beam. These time-series represent displacements and accelerations at a quarter of the beam and at the middle of the beam. It can be seen that a good agreement is obtained between measured and calculated signals of the uncontrolled piecewise linear beam system, i.e. (1) with u 0. Dierences can be seen more clearly in the acceleration signals where the shortcomings of the model are emphasized. A time-domain validation at various excitation frequencies is depicted in Fig. 3, which shows the maximum absolute values of the long-term periodic displacements at dierent excitation frequencies fexc . The good correspondence between numerical and experimental results illustrates the quality of the 3-DOF model. The long-term periodic behavior is also shown for the linear model, thus without the one-sided spring. There are two main dierences in periodic behavior between the linear and the nonlinear model. Firstly, the eigenfrequencies of the linear model precede the harmonic resonance frequencies of the nonlinear model; indicated with 1, due to the additional stiness of the one-sided spring in the nonlinear model. Secondly, coexisting n-periodic solutions n 2 f2; 3; . . .g with large vibration amplitudes appear for the nonlinear model. It can be seen that when solutions coexist, the 1-periodic solution is favorable with respect to vibration amplitude despite the fact that it sometimes becomes unstable; indicated with a dashed line. This justi®es the objective to stabilize 1-periodic solutions globally to achieve vibration amplitude reduction. To stabilize 1-periodic solutions of the uncontrolled beam system, naturally these solutions must exist within the frequency-range of interest. Fig. 3 shows that the existence of 1-periodic solutions can be shown
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66 –4
qact [m]
4
–4
quarter of the beam
x 10
4
numerical
qmid [m]
52
0 experimental
20 experimental numerical 0
–15 0
experimental
0 numerical
0.1
20 experimental numerical 0
–15 0
0.1
time [s]
middle of the beam
–12 0
0.1
q¨mid [ms–2]
q¨act [ms–2]
–12 0
x 10
0.1
time [s]
Fig. 2. 2-periodic solutions of (1) at fexc 35 Hz with u 0 N and qact at a quarter of the beam.
–2
–2
10 1
2
|qmid|max [m]
|qact|max [m]
10
3
–3
10
4 2
–4
10
o = experiments linear beam system
–5
10
1
10
10
2
3
–3
10
4
2 –4
10
o = experiments linear beam system
–5
2
fexc [Hz]
1
10
1
10
2
fexc [Hz]
10
Fig. 3. Periodic solutions of (1) at various excitation frequencies with u 0 N and qact at a quarter of the beam.
numerically or experimentally, the latter only in the case of stable solutions. However, it is also possible to prove the existence of 1-periodic solutions mathematically. In fact, the existence of 1-periodic solutions for the considered uncontrolled beam system can be derived for any uniformly bounded periodic excitation force applied anywhere between the middle of the beam and the leaf springs [10]. This results from the fact that solutions of (1) with u 0 can be proved to satisfy a uniform Lipschitz condition, i.e. there exists uniqueness of solutions. Furthermore, these solutions can proved to be ultimately bounded for any uniformly bounded excitation force [13]. Both uniqueness and ultimate boundedness of solutions are needed to obtain the mathematical proof for the existence of 1-periodic solutions of (1) with u 0, as is shown by Yoshizawa [24], see [10]. The existence of 1-periodic solutions is required for the study of the error dynamics related to these solutions. The error dynamics will be incorporated in the control design based on linear error feedback. 3. Linear error feedback Linear error feedback will be applied at a single location on the beam. So, the multi-DOF beam will be controlled with one control force only. Consequently, we are dealing with the so-called problem of underactuation [23]. Under-actuation, in general, prohibits a full linearization of the beam dynamics by feedback. Basically, this means that only part of the beam dynamics can be linearized by active control, whereas, in general, part of the beam dynamics will remain nonlinear. These so-called internal dynamics determines whether the entire system will behave as desired, i.e. whether the entire system will end up in the 1-periodic response. Linear error feedback does not aim at linearizing part of the beam dynamics. On the contrary, it
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66 2–periodic at 38 [Hz]
3–periodic at 57 [Hz] 200
φ(emid,t) [N]
φ(emid,t) [N]
300 knl•emid 0•emid 0
–100 –4
emid [m]
0
knl•emid 0•emid 0 –70 –3
2 x 10
0
1.5 –3
x 10
5
x 10
Φ(emid,t) [Nm–1]
Φ(emid,t) [Nm–1]
emid [m]
–3
5
2.5
53
knl
0
–0.5 0
1/38
time [s]
1/19
2.5
x 10
knl
0
–0.5 0
1/57
time [s]
1/19
Fig. 4. Sector-bounded nonlinearity /
emid ; t.
tries to change the beam dynamics as a whole by adding servo-damping or servo-stiness to the error dynamics of the system such that these dynamics, i.e. the underlying 1-periodic solutions, become locally, or globally, asymptotically stable. The error dynamics of the piecewise linear beam system are obtained by subtracting the equations of motion at the 1-periodic solution 1 p p
t, or Mp Bp_ Kp knl
pmid h2 hT2 p h2 v
4
T
with p pact pmid pn , from the equations of motion at an arbitrary solution q; see (1). This gives Me B_e Ke h2 /
emid ; t h1 u
5
T
with e q ÿ p eact emid en . A special feature of these error dynamics ± a feature which will be exploited in the stability analysis of these dynamics ± is the sector-bounded nonlinearity /
emid ; t, or /
emid ; t knl
emid pmid
t
emid pmid
t ÿ knl
pmid
tpmid
t:
6
The nonlinearity /
emid ; t is sector-bounded as it satis®es a so-called sector condition [15,21], or 0 6 U
emid ; t
/
emid ; t 6 knl ; emid
U
0; t 0:
7
Sector-boundedness results from the fact that the restoring force of the one-sided spring can be encapsulated by two linear spring forces: one with zero stiness, and one with stiness knl . The sector-boundedness of /
emid ; t is depicted in Fig. 4 for two periodic error solutions of (5) with u 0: a 2-periodic solution at 38 Hz, and a 3-periodic solution at 57 Hz; see also Fig. 3. It can be seen that /
emid ; t indeed satis®es the sector condition, i.e. U
emid ; t remains bounded by 0; knl . Note that U
emid ; t is periodic with a fundamental frequency xf 19 Hz whereas the excitation frequency equals to 38 and 57 Hz, respectively. So, U
emid ; t is also 2-periodic and 3-periodic, respectively. Obviously, the existence of periodic solutions for (5) shows that in the considered case no global stability for the ®xed point: e e_ 0, is obtained. Therefore, the error dynamics will be altered by choosing the input u in (5) as
1
Throughout this paper, desired long-term solutions will be approximated, if necessary, by means of truncated Fourier series [8].
54
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66
Fig. 5. Graphical representation of linear error feedback.
u ÿkd e_ act ÿ kp eact ;
8
with kp and kd positive ®nite control parameters. This input can be interpreted as a linear spring-damper combination with a prescribed motion at one end and the motion of the beam at the other end; see Fig. 5. (8) transforms (5) into Me Bc e_ Kc e h2 /
emid ; t 0
9
with Bc B kd h1 hT1
and
Kc K kp h1 hT1 :
10
Consequently, for all values of kp and kd , the stiness matrix Kc and the damping matrix Bc remain positive de®nite, hence Bc P B and Kc P K. In fact, similar to the uncontrolled case [10], all solutions of the error dynamics can be shown to be ultimately bounded irrespective of the values for the control parameters, i.e. solutions of the uncontrolled error dynamics cannot become unbounded by control. The error dynamics (9) may be locally, or globally, asymptotically stable depending on the values of the control parameters. The control parameters can be restricted to lower bounds and stability domains in the control parameter space. Lower bounds based on local asymptotic stability are obtained by the occurrence of bifurcations. Stability domains based on global asymptotic stability are derived with the circle criterion. 4. Local asymptotic stability The ability of linear error feedback to achieve local asymptotic stability will be investigated at dierent actuator locations and excitation frequencies for the ®xed point: e e_ 0. If the uncontrolled beam system, (5) with u 0, is excited at 11.7 Hz, the long-term behavior appears to be chaotic. This can be seen in the lower left part of Fig. 6, where one of the Lyapunov exponents converges to a positive value [16]. When the control parameters are gradually increased, the chaotic nature of the response will eventually be lost and, ®nally, local asymptotic stability of the 1-periodic response will be obtained. This can be seen in the bifurcation diagrams in the upper part of Fig. 6, where a reversed cascade of period-doubling bifurcations can be observed for increasing values of the control parameters kp and kd ; the bifurcation diagrams are obtained by monitoring the long-term velocity q_ act every period time of the excitation force. The chaotic response appears to be more sensitive to changes in kd compared to changes in kp . Hence increasing kd almost immediately results in stable n-periodic solutions, while increasing kp preserves the chaotic attractor in a large region. Local asymptotic stability of 1-periodic solutions can be related to the occurrence of a supercritical ¯ip bifurcation [18]. The ¯ip bifurcation determines a lower bound in the control parameter space; see the lower
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66 –3
–3
kd = 0
x 10
10
qact [ms–1]
qact [ms–1]
10
0
–5 0
0
15
[Nm–1]
30
42.1
50
kd [Nsm–1]
10
50 42.1
λ1≈+1.5
0
kd [Nsm–1]
Lyapunov exponents λi
kp
kp = 0
x 10
–5 0
5224 6000
55
λ2≈13.9 λ3≈17.2 λ4≈18.1 λ5≈20.1 λ6≈25.7
–30 0
0 0
170
5224 6000
kp [Nm–1]
time [s]
Fig. 6. Chaotic solutions of (9) at fexc 11.7 Hz with qact at a quarter of the beam.
right part of Fig. 6. Below this bound, 1-periodic solutions are unstable so other long-term solutions must coexist. If the uncontrolled beam system is excited at 35 Hz, the unstable 1-periodic solution of small vibration amplitude coexists with a stable 2-periodic solution of large vibration amplitude; see Fig. 3. During the process of increasing the control parameters, the amplitude of the 2-periodic error solutions decreases. At certain levels of the control parameters, 2-periodic solutions disappear due to the occurrence of a ¯ip bifurcation beyond which 1-periodic solutions become locally asymptotically stable. This is shown in Fig. 7, where the maximum absolute values of the displacements of 2-periodic error solutions are depicted for dierent values of the control parameters kp and kd .The occurrence of a ¯ip bifurcation can again be
–3
–3
|emid|max [m]
|eact|max [m]
x 10 1
0 0
0
x 10 1
0 0
0
4
kd [Nsm–1] 1000
4
x 10
6
–1 kp [Nm–1] kd [Nsm ] 1000
x 10 6
kp [Nm–1]
kd [Nsm–1]
|eξ|max [m]
600
0.01
0 0
0 4
kd [Nsm–1] 1000
x 10 6
kp [Nm–1]
0 0
kp [Nm–1]
5 4
x 10
Fig. 7. 2-periodic solutions of (9) at fexc 35 Hz with qact at a quarter of the beam.
56
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66 –4
–4
x 10
|emid|max [m]
|eact|max [m]
x 10 5
0
0
8
0
0
4
4
x 10
kd [Nsm–1]
500 6
x 10
–1 kp [Nm–1] kd [Nsm ]
kp [Nm–1]
600
0 0
= experimental
= experimental
= numerical
= numerical
kp [Nm–1]
kd [Nsm–1]
kd [Nsm–1]
600
500 6
0 0
5
kp [Nm–1]
4
x 10
5 4
x 10
Fig. 8. 2-periodic solutions experimentally determined at fexc 35 Hz with qact at a quarter of the beam.
regarded as a lower bound on the control parameters for the coexistence of solutions. This is shown in the lower right part of Fig. 7, where lines of equal error level are depicted together with values for the control parameters at the ¯ip bifurcation; these values are represented by + symbols. In general, the lower bound may be found using an ecient path following the method to track the bifurcation for varying control parameters [7]. The lower bound can also be determined experimentally as shown in Fig. 8. It can be seen that a broad resemblance with the behavior as shown in Fig. 7 results from the 20 20 experiments, i.e. the maximal absolute value of the displacements of 2-periodic error solutions decreases for increasing values of the control parameters until the lower bound, at which ¯ip bifurcations occur, is reached. Beyond the lower bound, 1-periodic solutions remain locally asymptotically stable. The occurrence of ¯ip bifurcations can also be used to assess the ability of linear error feedback to locally asymptotically stabilize 1-periodic solutions at dierent actuator locations along the beam. Hereto an equidistant grid of actuator locations is de®ned on the beam. According to this grid, the location at the left leaf spring is marked with 2 whereas the middle of the beam is marked with 67; see also Fig. 1. In between, 64 possible actuator locations are de®ned. At each location, the matrices M, B, and K dier because at each location the physical DOF qact is chosen dierently. In Fig. 9, the values of the control parameters at the ¯ip bifurcation are shown for various actuator locations: the left part for values of kp with kd 0 Nsmÿ1 , the right part for values of kd with kp 0 Nmÿ1 . It can be seen that 1-periodic solutions at an excitation kd = 0
5
kd [Nsm–1]
kp [Nm–1]
7
10
6
10
5
10
4
10
2
34
location of qact
67
10
kp = 0
4
10
3
10
2
10
2
34
location of qact
Fig. 9. Flip bifurcations of (9) at fexc 35 Hz for various actuator locations.
67
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66 location 34, kd = 0
0 25
|µi|max
fexc [Hz]
↔
100 1
2.5 4 0 25
fexc [Hz]
100 1
1
log(kp) [Nm–1]
8
log(kd) [Nsm–1]
|µi|max
8
2.5
↔
57
kp=2 106
1.2 1.4
1.8
1.6
2 2.2
1 25 30
56
70
fexc [Hz]
100
location 34, kp = 0
4
1 1.4 1.8 2
1 25 30
56
100
fexc [Hz]
Fig. 10. Local asymptotic stability of 1-periodic solutions of (9) with qact at location 34.
frequency of 35 Hz can be locally asymptotically stabilized at almost any actuator location ranging from the leaf springs, at location 2, to the middle of the beam, at location 67; note that the system is symmetric around the middle of the beam. However, away from the middle of the beam larger values for the control parameters are required due to the increased eective stiness of the beam. The ability to locally asymptotically stabilize 1-periodic solutions within a large frequency-range, and at a ®xed actuator location, is shown in Fig. 10. In this ®gure, linear error feedback is applied at location 34, which corresponds to a quarter of the beam; see Fig. 1. Various settings for the control parameters: kp in the upper part and kd in the lower part, are shown for a range of excitation frequencies. In the left part of Fig. 10, the local asymptotic stability of the 1-periodic solutions is shown represented by the maximum absolute Floquet multipliers jli jmax ; local asymptotic stability is obtained when jli jmax is less than one. In the right part of Fig. 10, lines of equal height with respect to jli jmax are shown. The lines denoted with 1 correspond to the transition from locally stable to locally unstable, or vice versa. This transition is represented by ¯ip bifurcations. It can be seen that starting at a frequency where 1-periodic solutions are unstable in the uncontrolled case ± roughly between 30 and 56 Hz ± ultimately gives local asymptotic stability if the control parameters are increased. Unfortunately, it can also be seen in the upper right part of Fig. 10, that at excitation frequencies where 1-periodic solutions initially are locally asymptotically stable, for example at 70 Hz, these solutions can become unstable for large values of the control parameter kp . This phenomenon also occurs for the control parameter kd , for example, for larger values of kd , closer to the leaf springs, or for smaller values of the modal damping f. Physically, linear error feedback results in the approximation of a node at the actuator location in the error dynamics for large values of the control parameters. This is shown in Fig. 11, where it can be seen that, although the actuator location approximately vibrates in the 1-periodic response, pact qact , the other DOFs vibrate in a 2-periodic response, pmid 6 qmid and pn 6 qn . Fig. 11 is based on an excitation frequency of 70 Hz whereas the control parameters are chosen as: kp 2 106 Nmÿ1 and kd 0 Nsmÿ1 , see the upper right part of Fig. 10. It can be seen that the control force does not converge to zero because the obtained long-term response is no natural response of the uncontrolled system. In fact, the internal force of the actuated DOF with the other DOFs requires control force to keep the actuator location vibrating close to the 1-periodic response. As this response, although hardly visible, is 2-periodic, the control force is also 2-periodic; note that u ÿkp eact .
58
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66
Fig. 11. Failure in stabilizing 1-periodic solutions of (9) at fexc 70 Hz with qact at location 34.
From the above, it may be concluded that the closed-loop system can still have coexisting solutions. Only if 1-periodic solutions would be globally asymptotically stable, the coexistence of other long-term solutions could be excluded. Therefore, global asymptotic stability of (9) will be investigated. 5. Global asymptotic stability based on the circle criterion A sucient condition for global asymptotic stability of (9) can be obtained with the circle criterion [13,15], which can also be related to stability domains in the control parameter space. Unfortunately, the circle criterion only guarantees global asymptotic stability for a limited number of actuator locations along the beam. To associate with the circle criterion, (9) is reformulated by the following state-space model _ A ÿ b2 w; y cT emid ;
11
w ÿ /
emid ; t with eT e_ T T , and 0 I ; A ÿM ÿ1 Kc ÿM ÿ1 Bc
h2 ; c 0
0 b2 : ÿM ÿ1 h2
12
This model constitutes the linear transfer between the input w and the output y, represented by the frequency response function v
jx: v
jx
Y
jx ÿ1 cT
A ÿ jxI b2 : W
jx
13
If A has no eigenvalues on the imaginary axis, if the nonlinearity /
emid ; t is sector-bounded, and if the linear error dynamics _ A ÿ b2 w; y cT emid ; w ÿ ki emid ;
14 ki 2 0; knl
are globally asymptotically stable, then the circle criterion can be applied to investigate whether 0, i.e. the ®xed point corresponding to 1-periodic solutions, is globally asymptotically stable. The circle criterion guarantees this stability when the inequality Rfv
jxg > ÿ
1 knl
8x 2 R
15
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66
59
Fig. 12. Graphical representation of the circle criterion with qact at location 54.
is satis®ed; Rf g denotes the real part of the corresponding argument. The circle criterion is based on the linear part of the error dynamics represented by the system matrix A. This linear part can be in¯uenced by the control parameters kp and kd ; see (10) and (12). The graphical representation of (15) requires v
jx to remain to the right of the vertical line through the point
ÿ1=knl ; 0 in the complex plane. This can be seen in Fig. 12, where v
jx is depicted for an actuator placed at location 54; this is closer to the middle of the beam than in the experimental setup where the actuator is placed at location 34. At location 54, global asymptotic stability of the error dynamics cannot be proved with the circle criterion in the uncontrolled case, kp kd 0, because v
jx does not remain to the right of the vertical line through the point
ÿ1=knl ; 0; see the left part of Fig. 12. However, it can be seen that after applying control with either kp 107 Nmÿ1 and kd 0 Nsmÿ1 , or kp 0 Nmÿ1 and kd 103 Nsmÿ1 , global asymptotic stability of the error dynamics is proved because v
jx now remains to the right of this line. Note that this stability result applies to all excitation frequencies. Stability domains in the control parameter space are determined by the smallest values of the control parameters that after being increased satisfy (15), i.e. all values that guarantee global asymptotic stability on the basis of the circle criterion. These domains are shown in Fig. 13, indicated by the hatched area, for an actuator placed at location 54, location 60, and location 66. Inside the stability domains, it will be assumed ± it has not been proved that (15) will, or will not, be satis®ed ± that 1-periodic solutions remain globally asymptotically stable as the essential part for stability in (15) hardly changes anymore. Outside these domains, global asymptotic stability cannot be guaranteed based on the circle criterion. It can be seen in Fig. 13, that away from the middle of the beam larger values for the control parameters are needed to guarantee global asymptotic stability. Furthermore, it can be seen in the left part of Fig. 13 that increasing the control parameters not necessarily implies that global asymptotic stability may still be guaranteed. For instance, when kp is increased whereas kd is kept at a constant value of 700 Nsmÿ1 , it is possible to escape from the stability domain such that global asymptotic stability is no longer guaranteed by the circle criterion. The ability to globally asymptotically stabilize 1-periodic solutions at an arbitrary location along the beam is shown in Fig. 14, where the maximum stiness of the one-sided spring for which the circle criterion guarantees global asymptotic stability is depicted as a function of the actuator location. The control location 54
location 60
kd [Nsm–1]
1400
location 66
700
700
700
0 0
kp [Nm–1]
4 6
x 10
0 0
kp [Nm–1]
2 6
x 10
0 0
kp [Nm–1]
Fig. 13. Stability domains for global asymptotic stability.
2 6
x 10
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66 kp = 1012, kd = 0
8
10
6
|min ℜ{χ(jω)}|–1
|min ℜ{χ(jω)}|
–1
60
knl ≈ 2 105 [Nm–1]
10
4
10
2
52
kp = 0, kd = 106
8
10
6
knl ≈ 2 105 [Nm–1]
10
4
10
67
2
52
location of qact
67
location of qact
Fig. 14. Global asymptotic stability of (9) at dierent actuator locations.
parameters are: kp 1012 Nmÿ1 and kd 0 Nsmÿ1 in the left part, and kp 0 Nmÿ1 and kd 106 Nsmÿ1 in the right part; both kp and kd are chosen extremely large to assure that the part of v
jx that determines the stability result is hardly in¯uenced for increasing the control parameters. It can be seen that global asymptotic stability with the current stiness of the one-sided spring can be guaranteed for both control parameters from the middle of the beam at location 67 to location 52. For actuator locations to the left of location 52, global asymptotic stability cannot be guaranteed with the circle criterion in its present form. However, the region of actuator locations for which global asymptotic stability may be guaranteed increases when the stiness of the one-sided spring, knl , is decreased. For example, global asymptotic stability can be guaranteed for all actuator locations when knl 6 104 Nmÿ1 . The circle criterion provides a sucient but not necessary condition for global asymptotic stability of 1-periodic solutions. Therefore, it induces conservatism into the control design. The amount of conservatism can be quanti®ed by studying the existence of long-term solutions for the error dynamics. 6. Existence of long-term solutions for the error dynamics The existence of long-term solutions for the error dynamics (9) can be restricted by the eigenfrequencies of both linear system matrices: A for qmid 6 0, and A0 for qmid > 0; A0 is related to A by changing Kc in (12), or 0 I 0
16 ; Kc0 Kc knl h2 hT2 : A ÿM ÿ1 Kc0 ÿM ÿ1 Bc This is motivated in Fig. 15, where the uncontrolled error dynamics, (5) with u 0, are depicted in terms of the frequency response function v
jx. If long-term solutions coexist, then (15) certainly is not satis®ed, i.e. every periodic solution of the error dynamics should satisfy for at least one frequency x P1 y 1 nÿ1 cn d
x ÿ nxf ;
17 Rfv
jxg 6 ÿ ; v
jx P1 w knl nÿ1 cn d
x ÿ nxf
imaginary axis
x 10
knl ≈ 2 105 [Nm–1]
–6
1
imaginary axis
–4
1
ω10,2 0 ωf
χ(jω) ω10,1
–3 –2
–1/k nl 0
real axis
2 –4
x 10
x 10
exploded view ω20,1=ω20,2
0 ω10,2
χ(jω) ω30,2 ωf –5 –4
ω30,1
–1/k nl 0
real axis
Fig. 15. Existence of periodic solutions of (1).
4 –5
x 10
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66
61
where cij denotes the jth Fourier coecient belonging to the periodic signal i. The individual frequency components are given by v
jnxf
Y
jnxf cy wn ; W
jnxf cn
n 2 f0; 1; 2; . . .g:
18
From (17), it is concluded that coexisting periodic error solutions must have a frequency component to the left of the vertical line through the point
ÿ1=knl ; 0. In Fig. 15, it can be seen that the coexisting 2-periodic solution at 38 Hz (+), and the 3-periodic solution at 57 Hz (), both based on xf 19 Hz, indeed have a frequency component that lies to the left of this line. Namely, the fundamental frequency component xf . In Fig. 15, also the values of the frequency response function v
jx at the eigenfrequencies of the system matrices A and A0 are depicted: xi0;1 for A and xi0;2 for A0 , i 2 f1; 2; 3g. It can be seen that the coexistence of long-term solutions for the uncontrolled piecewise linear beam system is somewhat conservatively restricted by the eigenfrequencies x10;1 and x10;2 . Moreover, x10;1 and x10;2 represent a lower and an upper bound, respectively, that, at least, should include one frequency component of a coexisting long-term solution. The lower bound ± in general de®ned by xi0;1 ± is based on the assumption Rfv
jxi0;1 g 0;
i 2 f1; 2; 3g:
19
This assumption approximately holds true for proportionally ± weakly ± damped systems such as the piecewise linear beam system. The lower bound represents a sucient restriction because the coexistence of long-term solutions by the circle criterion is restricted at Rfv
jxg ÿ1=knl instead of Rfv
jxg 0. The upper bound ± in general de®ned by xi0;2 ± is based on the same assumption, only now applied to the complementary system, ÿ1
v0
jx cT
A0 ÿ jxI b2
20
Rfv0
jxi0;2 g 0;
21
or i 2 f1; 2; 3g:
This assumption implies that Rfv
jx10;2 g ÿ1=knl , which can be derived from so-called loop transformations as shown in Fig. 16. The upper left part of Fig. 16 is the graphical representation of (11) and can be transformed to a system based on the frequency response function v0
jx and the sector-bounded non-
χ( j ω) χ( jω) χ( jω)
knl knl
φ() φ()
χ( j ω)
φ()
Fig. 16. Loop transformations.
φ()
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66 knl ≈ 2 105 [Nm–1]
–5
imaginary axis
1
x 10
–6
2
a2i +(1/knl)ai+b2i =0
ω10,2
χ(jω)
0 ω20,1,
ω20,2
ω30,2 ω30,1 –1 –1
0
–1/k nl
imaginary axis
62
knl ≈ 1 106 [Nm–1]
ω10,2
a2i +(1/knl)ai+b2i =0 χ(jω)
0 ω20,1, ω20,2 ω30,2
–4 –2
1
real axis
x 10
ω30,1
0
–1/k nl
2
real axis
–5
x 10
–6
x 10
Fig. 17. Solutions of Rfv0
jxi0;2 g 0.
linearity /0
emid ; t. From these loop transformations, it can be concluded that both frequency response functions v
jx and v0
jx are related by v0
jx
v
jx : 1 knl v
jx
22
Now, suppose there exists a solution v
jxi0;2 ai jbi , or v0
j xi0;2
ai
1 knl ai knl b2i jbi 2
2 2
1 knl ai knl bi
:
23
Then, (21) implies ai
1 knl ai knl b2i 0:
24
This equation represents a circle in the complex plane through the points
ÿ1=knl ; 0 and (0,0) as can be seen in Fig. 17. The points at which this circle crosses the frequency response function v
jx correspond to solutions of (21). In the left part of Fig. 17, it can be seen that both v
jx10;2 and v
jx30;2 correspond to such crossing points. Apparently, the underlying assumption, (21), approximately holds true in this case. In the right part of Fig. 17, the eect of increasing the stiness of the one-sided spring to knl 106 Nmÿ1 ± which is ®ve times the original value ± is shown. Now, it can be seen that the coexistence of long-term solutions is not only limited to frequency components in x10;1 ; x10;2 but also enables frequency components in x30;1 ; x30;2 . The upper bound ± de®ned by xi0;2 ± also represents a sucient restriction because, for the piecewise linear beam system, the point (ÿ1=knl ,0) represents the most negative real value that can be reached by v
jxi0;2 . Consequently, ai does not equal ÿ1=knl but is always somewhat larger; see Fig. 17. The existence of long-term solutions for the error dynamics can be used to illustrate the conservatism induced by the circle criterion. This is shown in Fig. 18, where bounds based on global asymptotic stability
log(kp) [Nm–1]
kp = 3.483 6
log(kd) [Nsm–1]
location 54, kd = 0
8
106
lb 1 1.4 1.7 2.2
ub 1 25 30
56
fexc [Hz]
100
location 54, kp = 0
4 lb
2.6
lb kd = 6.44 102 1 1.4 1.7
ub
2
1 25 30
56
fexc [Hz]
Fig. 18. Local asymptotic stability of (9) with qact at location 54.
100
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66
63
are depicted together with bounds based on local asymptotic stability. In this ®gure, lines of equal height are shown with respect to the largest Floquet multipliers li in absolute value at the ®xed point: e e_ 0, for an actuator placed at location 54. If jli jmax is larger than one, 1-periodic solutions are unstable and longterm solutions for the beam dynamics coexist. Fig. 18 shows that if a 1-periodic solution is unstable in the uncontrolled case: kp kd 0, roughly between 30 and 56 Hz, it can be locally asymptotically stabilized by increasing kp or kd . Global asymptotic stability can be guaranteed with the circle criterion, i.e. for all excitation frequencies when kp P 3:483 106 Nmÿ1 and kd 0 Nsmÿ1 , or when kp 0 Nmÿ1 and kd P 6:44 102 Nsmÿ1 ; see the left and right part of Fig. 18, respectively. The dashed lines in Fig. 18 represent the lower (lb) and upper (ub) bounds for the existence of 2-periodic solutions with fundamental frequencies of 2xf fexc , i.e. x10;1 6 xf 6 x10;2 , or 2x10;1 6 fexc 6 2x10;2 . These bounds, obtained with the eigenvalues of the system matrices: A and A0 , are both functions of the control parameters. Hence A and A0 depend on the control parameters. As 2-periodic solutions exist up to kp 106 Nmÿ1 and kd 0 Nsmÿ1 in the left part of Fig. 18, and up to kp 0 Nmÿ1 and kd 4 102 Nsmÿ1 in the right part of Fig. 18, global asymptotic stability at all excitation frequencies cannot occur below these values for the control parameters; note that below these values, however, global asymptotic stability may still be guaranteed at speci®c excitation frequencies. So, in these situations conservatism induced when applying the circle criterion is limited to 3:5 times the value for kp , or 1:5 times the value for kd , and might be even smaller. 7. Control simulations and experiments For the piecewise linear beam system, control simulations and experiments will be performed on the basis of local asymptotic stability only. This is because the ability of linear error feedback to guarantee local asymptotic stability is less constrained compared to the ability to guarantee global asymptotic stability with the circle criterion. In fact, in the control experiments an actuator is placed at a quarter of the beam, or location 34, where global asymptotic stability for all excitation frequencies cannot be guaranteed at all. Control simulations are shown in Fig. 19 at two dierent settings for the control parameters: in the left part only damping is added, kd 6 102 Nsmÿ1 , in the right part only stiness is added, kp
–4
emid [m]
–15 8
–4
0
0.2
0.2
1.5
0 –20 0
0.2
1.5
time [s]
kp = 3.5 104, kp = 0
–4
1.5
0.2
1.5
0.2
1.5
x 10
0
0
–0.01 0 30
u [N]
–0.01 0 30
8
–8 0 0.01
1.5
0
x 10
–15
eξ [m]
eξ [m]
5 0
1.5
x 10
–8 0 0.01
u [N]
–4
kp = 0, kd = 6 102
qact [m]
x 10
emid [m]
qact [m]
5 0
0 –20 0
0.2
1.5
time [s]
Fig. 19. Control simulations at fexc 35 Hz with qact at location 34.
64
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66 kp = 3.5 104, kd = 0
–4
qact [m]
5
x 10
0 experimental numerical
qmid [m]
–15 0 –4 x 10 5
0.2
1
0 experimental numerical
–15 0
0.2
1
30
u [N]
experimental
0 –20 0
numerical 0.2
1
time [s] Fig. 20. Control experiment at fexc 35 Hz with qact at location 34.
3:5 104 Nmÿ1 . No control is applied before t 0:2 s, such that the beam system vibrates in a long-term 2-periodic response of large vibration amplitude. At t 0:2 s, control is switched on after which all DOFs start vibrating according to the 1-periodic response of small vibration amplitude. As the beam tends to the 1-periodic response, the control force tends to zero. A similar behavior can be observed in real-life experiments performed on the lab-scale setup of the beam system. This is depicted in Fig. 20, where control is applied after t 0:2 s with kp 3:5 104 Nmÿ1 and kd 0 Nsmÿ1 . It can be seen that both qact and qmid approximately vibrate in the 1-periodic response at t 1 s whereas the control eort has become small. In between 0.2 and 1 s, the similarity in convergence of both numerical and experimental displacements, qact and qmid , can be seen clearly. The resemblance between numerical and experimental results could not be improved signi®cantly by repeating the experiment. This is not very surprising. Hence Fig. 8 already revealed a dierence between numerical and experimental values for the control parameters needed to create a lower bound for local asymptotic stability. Consequently, the resulting stiness in both the numerical and experimental error dynamics diers. Therefore, it is expected that the numerical and experimental error convergence diers as well for equal control parameters. The control performance appears to be very good on the lower bound where 2-periodic solutions disappear due to the occurrence of ¯ip bifurcations; see the upper part of Fig. 21. The control performance can be characterized by the control eort needed to stabilize 1-periodic solutions, or by the settling time needed to reach the 1-periodic response. The control eort needed to stabilize 1-periodic solutions is shown in the upper left part of Fig. 21 by means of numerical simulation for various values in the control parameter space. At those values where local asymptotic stability is obtained, the squared control force integrated over the time needed to reach the 1-periodic response, tf , is depicted. It can be seen that minimum control eort is obtained for the control parameters: kp 2:56 104 Nmÿ1 and kd 2:67 102 Nsmÿ1 . This setting can be found at the lower bound. Away from the lower bound, the amount of control eort usually increases. The settling time, tf , needed to reach the 1-periodic response shows a dierent behavior in the control parameter space as can be seen in the upper right part of Fig. 21. Increasing the values for kd with kp 0 Nmÿ1 seems to result in smaller settling time. However, less clear trends can be observed when increasing both values for the control parameters. Surprisingly, small settling times also occur near the
2
tf [s]
15
65
0
0
∫0
tf
u2(t) dt
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66
0
0
267
kd [Nsm–1] 1000
2.56 6
–4
[Nm–1]
0.2
–0.01 0
0.2
2.56 6
–4
1.5
0
0
267
kd [Nsm–1] 1000
emid [m]
0
0 0.01
eξ [m]
kp
x 10
u [N]
qact [m]
x 10
0 4
1.5
time [s]
8
x 10
kp
4
x 10
[Nm–1]
0 –8 0 30
0.2
1.5
0 –20 0
0.2
1.5
time [s]
Fig. 21. Control performance at fexc 35 Hz with qact at location 34.
control parameter setting where minimal control eort is obtained. Thus at the lower bound where the necessary values for the control parameters are applied; this is also observed at other excitation frequencies [9]. In the lower part of Fig. 21, the controlled behavior is shown with linear error feedback based on this control parameter setting. This shows the improvement that may be obtained by ®ne-tuning the control parameters. 8. Conclusions The 1-periodic solutions of the piecewise linear beam system may be stabilized both locally and globally to obtain vibration amplitude reduction with limited control eort. Local asymptotic stability can be related to the occurrence of bifurcations. These bifurcations can be used to obtain lower bounds in the control parameter space. Global asymptotic stability may be guaranteed with the circle criterion, and can be related to stability domains in the control parameter space. Global asymptotic stability on the basis of the circle criterion is only guaranteed for a limited range of actuator locations along the beam. This range can be extended by decreasing the one-sided spring constant. The ability to stabilize 1-periodic solutions on the basis of local asymptotic stability only is shown by various simulations and experiments. It has been shown that good results in terms of control eort and settling time can be obtained at the lower bound where just enough control force is applied needed to create a bifurcation. However, even local asymptotic stability cannot be obtained at all combinations of actuator locations and excitation frequencies. These combinations reveal the limitations of the present control design. Appendix A For an actuator placed at a quarter of the beam, location 34, the following matrices were used: 2
2:49171266 100
6 M 4 ÿ7:39353071 10ÿ1 ÿ1:00670592 10ÿ2
ÿ7:39353071 10ÿ1
ÿ1:00670592 10ÿ2
3
4:70207355 100
7 2:55145832 10ÿ2 5kg;
2:55145832 10ÿ2
2:78179792 10ÿ4
66
M.F. Heertjes, M.J.G. Van De Molengraft / Chaos, Solitons and Fractals 12 (2001) 49±66
2
9:13975994 101
6 B 4 ÿ5:46796875 101 2
ÿ5:46796875 101
ÿ3:12421552 10ÿ1
3
7:30340017 101
7 3:54360794 10ÿ1 5Nsmÿ1 ;
ÿ3:12421552 10ÿ1
3:54360794 10ÿ1
1:08926256 10ÿ2
3:82896231 105
ÿ2:59390789 105
9:78184925 10ÿ8
6 K 4 ÿ2:59390789 105 9:78184925 10ÿ8
3
2:10586254 105
7 ÿ2:65589609 10ÿ7 5Nmÿ1 ;
ÿ2:65589609 10ÿ7
6:10747501 101
knl 1:966 105 Nmÿ1 ; c 9:398 10ÿ4 Ns2 :
References [1] Berg J, Briggs GAD. Nonlinear dynamics of intermittent-contact mode atomic force microscopy. Phys Rev B 1997;55(22):14899± 908. [2] Bishop SR. Philosophical transactions of the royal society london. Impact Oscillators 1994;347:347±51. [3] Chen G, Dong X. From chaos to order: perspectives and methodologies in controlling chaotic nonlinear dynamical systems. Int J Bifurc Chaos 1993;3:1363±409. [4] Craig RR, in: Martinez DR, Miller AK, editors. Combined experimental/analytical modeling of dynamic structural systems 1±30, New York: ASME Appl. Mech. AMD-67, A review of time-domain and frequency-domain component mode synthesis method, 1985. [5] Doole SH, Hogan SJ. A piecewise linear suspension bridge model: nonlinear dynamics and orbit continuation. Dyn Stab Syst 1996;11:19±47. [6] dSPACE digital signal processing and control engineering GmbH, User's guide, dSPACE GmbH, 1994. [7] Fey RHB, van Campen DH, de Kraker A. Long term structural dynamics of mechanical systems with local nonlinearities. ASME J. Vib Acous 1996;118:147±53. [8] Heertjes MF, van de Molengraft MJG, Kok JJ, van Campen DH. Vibration reduction of a harmonically excited beam with onesided spring using sliding computed torque control. Dyn Cont 1997;7:361±75. [9] Heertjes MF, van de Molengraft MJG. Enschede NOLCOS'98: IFAC Global stabilization of the 1-periodic response of a beam with one-sided spring. Fourth Nonlinear Control Systems Design Symposium 1998:47±52. [10] Heertjes MF, van de Molengraft MJG, Kok JJ. Proceedings of the 1999 ASME Design Engineering Technical Conferences, VIB8031 (1999) 1±7, Las Vegas, DETC99 Uniformly bounded input gives ultimately bounded output ± On the dynamics of a piecewise linear beam system. [11] Hu HY. Controlling chaos of a periodically forced nonsmooth mechanical system. Acta Mech Sinica 1995;11:251±8. [12] Karagiannis K, Pfeier F. Theoretical and experimental investigations of gear-rattling. Nonlinear Dynamics 1991;2:367±87. [13] Khalil HK, Nonlinear Systems, second edition, New Jersey: Prentice Hall, 1996. [14] Lazer AC, McKenna PJ. Large-amplitude periodic oscillations in suspension bridges: some new connections with nonlinear analysis. SIAM Rev 1990;32(4):537±78. [15] Leonov GA, Ponomarenko DV, Smirnova VB. Frequency domain methods for nonlinear analysis: theory and application. Singapore: World Scienti®c; 1996. [16] Parker TS, Chua LO. Practical numerical algorithms for chaotic systems. New York: Springer; 1989. [17] Pfeier F, Kunert A. Rattling models from deterministic to stochastic processes. Nonlinear Dynamics 1990;1:63±74. [18] Shaw SW, Holmes PJ. A periodically forced piecewise linear oscillator. J. Sound and Vibration 1983;90:129±55. [19] Shinbrot T. Progress in the control of chaos. Adv Phys 1995;44:73±111. [20] SIMULINK dynamic system simulation software, User's guide, The MathWorks Inc., 1995. [21] Slotine J-JE, Li W. Applied nonlinear control. New York: Prentice Hall; 1991. [22] Spatz JP, Sheiko S, M oller M, Winkler RG, Reineker P, Marti O. Forces aecting the substrate in resonant tapping force microscopy. Nanotechnology 1995;6:40±4. [23] Spong MW. The swing up control problem for the acrobot. IEEE Cont Syst Mag 1995:49±55. [24] Yoshizawa T. Stability theory and the existence of periodic solutions and almost periodic solutions. New York: Springer; 1975.