# Numerical Exploratory Analysis of Dynamics and Control of an Atomic Force Microscopy in Tapping Mode with Fractional Order

Authors: Mauricio A. Ribeiro; Jose M. Balthazar; Wagner B. Lenz; Rodrigo T. Rocha; Angelo M. Tusset
Journal: Shock and Vibration
 (2020)

---

## Abstract

In this paper, we investigate the mechanism of atomic force microscopy in tapping mode (AFM-TM) under the Casimir and van der Waals (VdW) forces. The dynamic behavior of the system is analyzed through a nonlinear dimensionless mathematical model. Numerical tools as Poincaré maps, Lyapunov exponents, and bifurcation diagrams are accounted for the analysis of the system. With that, the regions in which the system presents chaotic and periodic behaviors are obtained and investigated. Moreover, the fractional calculus is introduced into the mathematical model, employing the Riemann-Liouville kernel discretization in the viscoelastic term of the system. The 0-1 test is implemented to analyze the new dynamics of the system, allowing the identification of the chaotic and periodic regimes of the AFM system. The dynamic results of the conventional (integer derivative) and fractional models reveal the need for the application of control techniques such as Optimum Linear Feedback Control (OLFC), State-Dependent Riccati Equations (SDRE) by using feedback control, and the Time-Delayed Feedback Control. The results of the control techniques are efficient with and without the fractional-order derivative.


---

## Body

## 1. Introduction

Technological advances in the development of electromechanical systems are gaining ground in the most diverse branches of engineering science. Such advances permit the development of smaller devices that vary from macro- to nanoscale, which have opened space for new research fields. However, the size scale of these devices has been a challenge as classical mechanics is not the only applicable one anymore. Mainly, nanoelectromechanical systems (NEMS) are affected by quantum forces, whose systems have been extensively studied in the past years [1–4].A special mechanism has been used for sample surfaces analysis at atomic scale, which is the atomic force microscope characterized as a NEMS, mostly referred to as atomic force microscopy technique. This technique is very well established as it is a very precise superficial analysis and has allowed the increase of the understanding and analysis of very small and soft materials such as polymeric materials [5], ceramics [6], biological cells [7], and surface tribological analyses [8]. The mechanisms are found to be actuated by a piezoelectric material which is also used for controlling vibrations [9, 10].Among the AFM in tapping mode (AFM-TM), contact and noncontact with the sample surface stand out, as they can form a three-dimensional image of such surface [11]. According to [12], it has been experimentally observed that the microcantilever of AFM in tapping mode may undergo chaotic behavior under certain conditions, due to the interactions of atomic scale forces between the surface and the tip of the AFM-TM [12].The van der Waals (VdW) forces are the most dominant atomic force in the AFM-TM and can be found in various works in the literature as in [12]. The authors provide the mathematical modeling of the dynamical equations of motion of the AFM-TM system to investigate the chaotic behavior of the microcantilever beam of the AFM-TM mechanism under VdW forces, where chaos is found and a control design has to be proposed.The VdW forces consider the electron density fluctuations present between the test tip atoms and the surface. In [13], the authors proposed an AFM-TM model with hydrodynamic and damping forces for a thin film. In this case, the mathematical model is a dumped mass-spring-damped system in which the acting force is derived from the Lennard-Jones potential, where a chaotic process is encountered in the system [14, 15].On the other hand, different quantic forces can become as dominant as VdW depending on the material of the surface to be scanned by the AFM mechanism. The authors in [16] investigate the transition of this other force with the VdW one through the AFM-TR when using a macroscopic gold surface on a flat geometry, whose force is the Casimir force. The transition between these forces appears with 10% separations of the λp plasma wavelength for the gold surfaces. In this analysis, the authors estimate the Hamaker constant for the VdW force regime that is in agreement with the precision of the Lifshitz theory, considering the surface roughness corrections. The experiment was carried out in a Pico Force system, with a gold sphere of approximately 100 μm of diameter over a microcantilever of length l = 100 nm.The intersection of the transition between the Casimir and VdW forces is also discussed in [17]. The authors discuss when the Casimir force becomes as dominant as VdW forces depending on the distance between the AFM microcantilever and the sample surface, concluding that there is a coexistence of the forces at very short distances. Recently, an analytical estimation of the Casimir force is found in [18], where the authors use a circular microplate in a two-sided NEMS capacitive system. The authors modeled the NEMS system and analyzed it using the harmonic balance method, where stability conditions, bifurcation points, and frequency responses are accounted.The atomic forces have been under study and are of great interest due to their effect on the dynamics of the microcantilever beam of the AFM. The mechanism becomes inoperable due to irregular measurements as mentioned in [12–14]. In addition, in some circumstances, the Casimir force becomes very expressive and can be dominant along with the VdW forces. This interaction between the forces brings a great need for its understanding and study, as the technology fastly advances to smaller and smarter devices, where the devices’ size scale makes these forces dominant.In this work, we approach computationally the dynamics and control of the AFM-TM model proposed by [14], accounting for the addition of the Casimir force along with the van der Waals forces for the interaction of the tip with the surface of the sample. The addition of the Casimir force states implications on the resulting dynamics of the system.There is an analysis of the nanosystem for a particular case where there is the coexistence of both Casimir and VdW forces in the system. Since the Casimir force is a problem to be solved yet, it is considered as a function and only numerical simulations of the system are carried out. Consolidated numerical techniques as Poincaré maps, Lyapunov exponent calculus, and bifurcation diagram are carried out for the system [19, 20].In addition, squeeze-film damping is also introduced due to the small distance between the microcantilever and surface as a viscoelastic term. The viscoelastic term is an approximation for the behavior of the indented analysis tip during the tapping process, commonly observed in the analysis of biological samples that, due to high vibrations and small distances, a gas film is generated [21]. With that, a fractional-order derivative is introduced. The fractional model allows the analysis of the viscoelastic behavior of the system on the action of the aforementioned forces in a computational way. However, due to the difficulty in obtaining the Jacobian algorithms, to the system in fractional order, the 0-1 test is utilized [14].In both fractional model and nonfractional models, the chaotic behavior of the system is analyzed and control techniques are proposed in order to control the chaotic behavior in the dynamics of the AFM system, which are Optimum Linear Feedback Control (OLFC) and State-Dependent Riccati Equations (SDRE) feedback controls and the Time-Delayed Feedback Control (TDFC).The rest of the paper is organized as follows. Section2 describes the mathematical model and the forces acting on the system. Section 3 describes the proposal of the controllers, their design, and application in the system without the fractional-order derivative. Section 4 discusses the fractional model of the AFM-TM system and its dynamic behavior with relevant numerical techniques. Section 5 describes the control techniques applied for the fractional model. And, finally, the conclusions of the work are stated in Section 6.


## 2. Mathematical Model for a Typical AFM-TM with a Nonlinear Behavior

Figures1(a) and 1(b) show the schematics for an AFM-TM on tapping mode. It is represented as a mass-spring-damper system, excited by a harmonic force induced by a piezoelectric actuator. The interaction force between the tip of the cantilever and the surface of the sample can be simplified as the interaction of the needle sphere and the surface, represented by(1)Ux,z0=A1R1260z0+x7−A2R6z0+x,where z0 is the distance between the equilibrium point of the microcantilever and the analyzed sample, R is the radius of the sphere of the tip of the AFM-TM, U (x, z0) is the Lennard-Jones potential, A1 is the Hamaker constant to the attractive potential, and A2 is the Hamaker constant to the repulsive potential.Figure 1
Representation of AFM-TM schematic diagram. (a) Representation of microcantilever beam and (b) representation of analysis tip through a mass-spring-damper model.
(a)(b)The van der Waals forces is described as a combination of the attractive and the repulsive parcels yielding(2)FVdW=−∂U∂z0+x=A1R180z0+x8−A2R6z0+x2.As the AFM-TM analysis is carried out on tapping mode, there is an excitation force on the microcantilever by the piezoelectric actuator given byf0cos (ωt), guaranteeing the oscillatory contact generated between the tip and the surface of the sample. Consequently, the equation of motion of the microcantilever beam is given by(3)mx¨=Fk+Fc+Fcs+FVdW+Fcas+f0cosωt,where x is the displacement of the tip of the microcantilever, m is the mass of the microcantilever, FVdW is the attraction and repulsive forces described by (2), the spring force is Fk, the structural damping force is Fc, the squeeze-film damping is given by Fcs, and the Casimir force is Fcas.The conservative force of the springFk is given by(4)Fk=klx+knlx3,where kl is the linear stiffness and knl is the nonlinear stiffness. The structure damping force Fcis described by(5)Fc=cdx˙,where cd is the structural damping coefficient of the microcantilever. The damping force generated by the squeeze-film damping Fcs is denoted by(6)Fcs=μeffη3lz0+x3x˙,where μeff is the coefficient of effective viscosity and η and l are the width and length of the microcantilever, respectively.The Casimir forceFcas is given by(7)Fcas=hclπ2240z0+x4,where h is Planck’s constant and cl is the speed of light.Based on the forces described in (2)–(7), it is possible to rewrite (3) as(8)mx¨+cdx˙+klx+knlx3=A1R180z0+x8−A2R6z0+x2−μeffη3lz0+x3x˙+hclπ2240z0+x4+f0cosωt.Carrying out a dimensionless procedure into (8), the equations of motion of the microcantilever beam are rewritten in the nondimensional form as(9)x˙1=x2,x˙2=−rx2−bx1−cx13+ea+x18−da+x12−pa+x13x2+βa+x14+ΓcosΩτ,where the dimensionless coefficients are(10)x1=xzs,x2=x˙ω1zs,τ=ω1t,a=z0zs,b=klkat,c=knlkatzs2,d=427,b=2a6405zs6,Γ=f0katzs,p=μeffη6lmω1zs3,Ω=ωω1,zs=32A2R3kl1/3,ψ=hπ2cl240,r=1Q,ω1=km,β=ψω12zs5,where ω1 is the frequency of the first mode of vibration of the microcantilever and Q is the quality factor coefficient that depends on the viscous fluid medium where the AFM cantilever exists [13]. In addition, zsis the equilibrium distance variable and is defined by zs=3/22D1/3, where D=A2R/6kl.Table1 shows the values of the physical constants that are accounted for the numerical simulations. These values are approximate values of those used during the experimental analysis that were obtained by the authors in [13, 16]. The choice of these parameters is because they were obtained when there is an interaction between the Casimir and VdW forces.Table 1
Parameters of the AFM-TM model [13, 16].
DescriptionValueLength of the cantilever449μmWidth of the cantilever46μmThickness of the cantilever1.7 μmTip radius0.15μmMaterial density2,330 kg/m3Young’s modulus176 GPaBending stiffness0.11 Nm−1Frequency of the 1st mode of the cantilever11804 kHzQuality factor100Hamaker constant (repulsive)1.3596 × 10−70 J⋅m6Hamaker constant (attractive)1.865 × 10−19 JUsing Table1 into the dimensionless coefficients of (9) yields r = 0.1, a = 1.6, b = 0.05, c = 0.35, d = 4/27, e = 0.0001, Γ = 0.2, p = 0.009, and β ∈ [−0.1 : 0.25]. The parameters β and p are accounted for analysis as they describe the intensity of the Casimir force and the viscoelastic term of the system. The system is assumed to be excited in resonance with the 1st mode of vibration (Ω = 1.0) [13].


### 2.1. Influence of the Squeeze-Film Damping and Casimir Force in the Behavior of the AFM-TM System

In this subsection, the dynamic behavior influenced by the squeeze-film damping and the Casimir force is investigated. Numerical simulations are carried out by using the 4th order Runge-Kutta implicit method with integration steph = 0.01 (ode45 of Matlabr) [22].The squeeze-film damping is inherent to the sample, providing a nonlinear damping. However, the Casimir force is only presented at nanodistances [23]. In this case, the interaction sphere plane between the surfaces is repulsive.Figure2 depicts the highest Lyapunov exponent with a variation on the parameter p ∈ [0.001 : 0.009] and β ∈ [−0.2 : 0.2], where the Jacobian algorithm for the Lyapunov calculus is carried out as in [23]. The initial conditions are (0, 0) and a computational time for the convergence is t = 10000s with a transient time of approximately 40% of the convergence time, where positive exponent values are denoted from yellow to green color [0 : 0.15] and negative values are light gray to black [−0.3 : 0].Figure 2
Lyapunov exponent in 2D withr = 0.1, a = 1.6, b = 0.05, c = 0.35, d = 4/27, e = 0.0001, Γ = 0.2, p ∈ [0.001 : 0.009], Ω = 1.0, and β ∈ [−0.1 : 0.25].It is observed that, in the region of the parameters (β vs p), Lyapunov exponents are positive and show an evidence of chaotic behavior in the interval of β ∈ [0.001 : 0.009] for any value of p. Then, for further numerical simulations and discussions, p = 0.009 is adopted for the viscoelastic term.Figures3(a) and 3(b) show the bifurcation diagram and Lyapunov exponent, respectively, for the variation of β in the interval [−0.1 : 0.25], with the adopted p = 0.009. The ranges of β that yield positive high values of the higher Lyapunov exponent and present chaotic behavior are when β ≈ [−0.0027 : 0.0025], β ≈ [0.0095 : 0.0675], β ≈ [0.0677 : 0.0685], and β ≈ [0.0924 : 0.0926]. A local value of β = 0.01818 is chosen for further analyses due to the presence of chaotic behavior, highest positive Lyapunov exponent, and it is in a very sensitive region of transition of behavior from chaos to periodic behavior.Figure 3
(a) Bifurcation diagram and (b) Lyapunov exponent withr = 0.1, a = 1.6, b = 0.05, c = 0.35, d = 4/27, e = 0.0001, Γ = 0.2, p = 0.009, Ω = 1.0, and β ∈ [−0.1 : 0.25].
(a)(b)Figures4(a) and 4(b) present the Poincaré map and the phase plane, respectively, for β = 0.01818 and p = 0.009. The existence of the chaotic behavior of the system is clear.Figure 4
(a) Poincaré map and (b) phase planet withr = 0.1, a = 1.6, b = 0.05, c = 0.35, d = 4/27, e = 0.0001, Γ = 0.2, p = 0.009, Ω = 1.0, and β = 0.01818.
(a)(b)For a temporal view of the chaotic behavior forβ = 0.01818 and p = 0.009, Figures 5(a) and 5(b) depict the time histories of displacement and velocity of the microcantilever beam, respectively, considering the time interval of 600≤τ≤1000 after transient period.Figure 5
(a) Time history of displacementx1 and (b) time history of velocity x2 with r = 0.1, a = 1.6, b = 0.05, c = 0.35, d = 4/27, e = 0.0001, Γ = 0.2, p = 0.009, Ω = 1.0, and β = 0.01818.
(a)(b)With the previous analyses of the system that showed the presence of chaotic behavior for parametersr = 0.1, a = 1.6, b = 0.05, c = 0.35, d = 4/27, e = 0.0001, Γ = 0.2, p = 0.009, Ω = 1.0, and β = 0.01818, the next section presents the design of different control techniques in order to suppress and lead the system to a customized desired periodic orbit, which is of great interest for the AFM application.


## 2.1. Influence of the Squeeze-Film Damping and Casimir Force in the Behavior of the AFM-TM System

In this subsection, the dynamic behavior influenced by the squeeze-film damping and the Casimir force is investigated. Numerical simulations are carried out by using the 4th order Runge-Kutta implicit method with integration steph = 0.01 (ode45 of Matlabr) [22].The squeeze-film damping is inherent to the sample, providing a nonlinear damping. However, the Casimir force is only presented at nanodistances [23]. In this case, the interaction sphere plane between the surfaces is repulsive.Figure2 depicts the highest Lyapunov exponent with a variation on the parameter p ∈ [0.001 : 0.009] and β ∈ [−0.2 : 0.2], where the Jacobian algorithm for the Lyapunov calculus is carried out as in [23]. The initial conditions are (0, 0) and a computational time for the convergence is t = 10000s with a transient time of approximately 40% of the convergence time, where positive exponent values are denoted from yellow to green color [0 : 0.15] and negative values are light gray to black [−0.3 : 0].Figure 2
Lyapunov exponent in 2D withr = 0.1, a = 1.6, b = 0.05, c = 0.35, d = 4/27, e = 0.0001, Γ = 0.2, p ∈ [0.001 : 0.009], Ω = 1.0, and β ∈ [−0.1 : 0.25].It is observed that, in the region of the parameters (β vs p), Lyapunov exponents are positive and show an evidence of chaotic behavior in the interval of β ∈ [0.001 : 0.009] for any value of p. Then, for further numerical simulations and discussions, p = 0.009 is adopted for the viscoelastic term.Figures3(a) and 3(b) show the bifurcation diagram and Lyapunov exponent, respectively, for the variation of β in the interval [−0.1 : 0.25], with the adopted p = 0.009. The ranges of β that yield positive high values of the higher Lyapunov exponent and present chaotic behavior are when β ≈ [−0.0027 : 0.0025], β ≈ [0.0095 : 0.0675], β ≈ [0.0677 : 0.0685], and β ≈ [0.0924 : 0.0926]. A local value of β = 0.01818 is chosen for further analyses due to the presence of chaotic behavior, highest positive Lyapunov exponent, and it is in a very sensitive region of transition of behavior from chaos to periodic behavior.Figure 3
(a) Bifurcation diagram and (b) Lyapunov exponent withr = 0.1, a = 1.6, b = 0.05, c = 0.35, d = 4/27, e = 0.0001, Γ = 0.2, p = 0.009, Ω = 1.0, and β ∈ [−0.1 : 0.25].
(a)(b)Figures4(a) and 4(b) present the Poincaré map and the phase plane, respectively, for β = 0.01818 and p = 0.009. The existence of the chaotic behavior of the system is clear.Figure 4
(a) Poincaré map and (b) phase planet withr = 0.1, a = 1.6, b = 0.05, c = 0.35, d = 4/27, e = 0.0001, Γ = 0.2, p = 0.009, Ω = 1.0, and β = 0.01818.
(a)(b)For a temporal view of the chaotic behavior forβ = 0.01818 and p = 0.009, Figures 5(a) and 5(b) depict the time histories of displacement and velocity of the microcantilever beam, respectively, considering the time interval of 600≤τ≤1000 after transient period.Figure 5
(a) Time history of displacementx1 and (b) time history of velocity x2 with r = 0.1, a = 1.6, b = 0.05, c = 0.35, d = 4/27, e = 0.0001, Γ = 0.2, p = 0.009, Ω = 1.0, and β = 0.01818.
(a)(b)With the previous analyses of the system that showed the presence of chaotic behavior for parametersr = 0.1, a = 1.6, b = 0.05, c = 0.35, d = 4/27, e = 0.0001, Γ = 0.2, p = 0.009, Ω = 1.0, and β = 0.01818, the next section presents the design of different control techniques in order to suppress and lead the system to a customized desired periodic orbit, which is of great interest for the AFM application.


## 3. Design of the Proposed Controls

### 3.1. Time-Delayed Feedback Control (TDFC)

As originally suggested by the author in [25], a continuous control input to stabilize a chaotic oscillation is given by the difference between the current output and the past one as follows [12, 25–27]:(11)U=κgx1τ−T,x2τ−T−gx1τ,x2τ,where T is the time delay and κ is the feedback gain. The term gx1τ,x2τimplies a scalar output signal measured at the current time T and the past time τ−T, respectively. Since the control input equation (11) only depends on the output signal, the time delay T is adjusted to the period of a targeted stable periodic orbit that is intended to be stabilized in a chaotic attractor. The control input, therefore, converges to null after the controlled system is stabilized to the targeted orbit.Assuming that the velocity of oscillation is measured as an output of the nonlinear system (9), the control signal is given by(12)U=κx1τ−T,x2τ−T.The NEMS system with the control signal of (12) is expressed in the following way:(13)x˙1=x2,x˙2=−rx2−bx1−cx13+ea+x18−da+x12−pa+x13x2+βa+x14+ΓcosΩτ+U.The time delayT and feedback gain κ are important control parameters that substantially affect the control performance. The time delay T is adjusted to T=2π/Ω=2π/1=2π to stabilize an orbit with the same frequency as the external force oscillating the system. Figure 6 shows a bifurcation diagram of the system related to the gain κ.Figure 6
Bifurcation diagram for control gainκ.For anyκ≥0.06, the controlled system has a periodic behavior, with unitary period. Then, the gain is defined as κ=1.7 which is a value far from the chaotic region and that guarantees that the system will remain in the periodic orbit with unitary period. Figure 7 shows the comparison of the uncontrolled system (9) to the controlled one (13).Figure 7
Time histories of (a) displacement and (b) velocity; (c) phase portrait of the system with controlled trajectory (solid line) and the trajectory without control (dashed line); (d) control signalU.
(a)(b)(c)(d)The obtained periodic orbits with the time delay control of Figure7 can be numerically calculated through a Fourier series as(14)x1=−0.61128654−0.26404525sinτ,x2=−0.26404525cosτ.The time-delayed control was efficient in leading the system to a periodic orbit. As can be seen in the results, the control only used the control signal when necessary to take the system to a periodic orbit of the system, as observed in Figure7(d), where the control signal (U) tends to zero as the system stabilizes in the periodic orbit (Figures 7(a) and 7(b)). These results can also be observed in [12, 14].However, for the cases where it is desired to take the system to a previously defined orbit, the control may not be the most suitable. For these mentioned cases, the optimal control is considered in the next section.


### 3.2. Optimal Control

In the following, the solution of these problems using the optimal linear control technique for nonlinear systems developed by the authors in [28] is presented. The introduction of the control signal in (9) for the optimal control is considered equally in (13).However, the vector controlU for this optimal control consists of two parts: U=u˜+u, where u is the linear feedback control and u˜ is the feedforward control, and the control that maintains the system in the desired trajectory is given by(15)u˜=−x˜⋅2+rx˜2+bx˜1+cx˜13−ea+x˜18+da+x˜12+pa+x˜13x˜2−βa+x˜14−ΓcosΩτ.Substituting (15) into (13) and defining the deviation of the desired trajectory, (16)e=x1−x˜1x2−x˜2,where x˜1τand x˜2τ are the trajectories defined by (14).The system can be represented in the following form:(17)e˙1=e2,e˙2=−re2−be1−ce1+x˜13+ea+e1+x˜18−da+e1+x˜12−pa+e1+x˜13x˜2+βa+e1+x˜14+cx˜13−ea+x˜18+da+x˜12+pa+x˜13x˜2−βa+x˜14+u.


#### 3.2.1. Feedback Control by OLFC

The optimal linear feedback control is applied as it was introduced by the authors in [28]. The system of (17) represented in deviations form can be written in matrix form as(18)e˙=Ae+Ge,x˜+Bu.According to [12], if there are matrices Q and R, positive and definite, with Q symmetric matrix, such that(19)Q˜=Q−GTe,x˜P−PGe,x˜,being positive and definite, in which the matrix G is limited, then the control u is optimal and it transfers the nonlinear system of (18) from any initial state to the final state e∞=0.Minimizing the cost functional to(20)J=∫0∞eTQ˜e+uTRudτ.The controlu can be found solving the following equation:(21)u=−R−1BTPe.HavingP as a symmetric matrix, the algebraic Riccati equation is developed, denoted by(22)PA+ATP−PBR−1BTP+Q=0.The matricesA and B may be represented by(23)A=01−b−r,B=01,and by definition(24)Q=10000100,R=1.Solving (22), the values of the cells of the matrix P are obtained, given by(25)P=109.49999.9019.90110.8503.Substituting the matricesR, B, and P into (21), the OLFC control signal, in deviation form, is given by(26)u=−9.9501e1−10.85504e2=−9.9501x1−x˜1−10.85504x2−x˜2.In Figure8, we observed the controlled system of (13) with control (U=u˜+u) in the designed orbit ((14).Figure 8
(a) Displacement of the tip with control (solid line) and without control (dashed line); (b) phase portrait of the tip with control (solid line) and without control (dashed line); (c) error variatione1=x1−x˜1 and e2=x2−x˜2.
(a)(b)(c)According to Figure8(c), the optimal linear feedback control was demonstrated to be effective in leading and maintaining the system in the desired orbit. Comparing the results obtained with the TDFC and OLFC controls, the OLFC control leads the system to the desired orbit much faster than the TDFC control. This is one of the advantages of using state feedback control (such as OLFC and SDRE). However, the determination of a desired orbit is needed, and, to keep the system in the desired orbit, the control signal has to be kept always acting in the system [12], which can be a disadvantage depending on the application of the controller. Nevertheless, a more suitable periodic orbit could be chosen depending on the purpose of the applied control.


#### 3.2.2. Feedback Control by SDRE

For the SDRE control, (17) may be represented in matrix form as(27)e˙=Aee+Ge,x˜+Bu,where the matrix of states A(e) also depends on the deviations e.Considering that, the matricesA and B may be represented by(28)A=01−b−ce12+3e1x˜1+3x˜12−r−p/a+e1+x˜13,B=01.In addition, by definition(29)Q=10000100,R=1.The quadratic performance measured for the feedback control problem is given by(30)J=12∫t0∞eTQe+uTRudτ,where Q and R are positive definite matrices. Assuming full state feedback, the control law is given by(31)uf=−R−1BTPee=−kee.HavingPe as a symmetric matrix, the algebraic Riccati equation is denoted by(32)PeAe+AeTPe−PeQR−1BTPe+Q=0,keis obtained for each iteration by solving Riccati equation (32), and, considering the asymptotically stable periodic orbit in (14) as desired orbits (xj=xj∗), the optimal (or suboptimal) control problem (u) can be formulated as follows: determine the control signal uf ((31)) that transfers (9) from the initial state(33)x10=x10,x20=x20,into the final state(34)x1∞=x˜1,x2∞=x˜2,and minimizes the functional (30). Such minimization shown in (30) implies the minimization of the system deviation ((27)) of the desired state (e=x−x˜) and of the state of the control applied to the feedback control (u).Another important factor to consider is that the matrixA(e) cannot violate the controllability of the system. The system of (13) is controllable if the rank of the matrix M is 2, according to the rule(35)M=B2×1A2×2eB2×1.Then, to obtain a suboptimal solution for the dynamic control problem, the SDRE technique has the following procedure:(1)
Define the state-space model with the state-dependent coefficient as in (27).(2)
Definex(0) = x0, so that the rank of M is n and choose the coefficients of weight matrices Q and R.(3)
Solve Riccati equation (32) for the state xτ.(4)
Calculate the input signal from (31).(5)
Integrate system (27) and update the state of the system xτ with these results.(6)
Calculate the rank of (35); if rank = 2, go to step 3. However, if rank <2, the matrix A(e) is not controllable; therefore, the last matrix controllable A(e) that has been obtained should be used, and thus go to step 3.Figure9 shows system (13) without control and with control (U=u˜+u) in designed orbit equation (14).Figure 9
(a) Displacement of the tip with control (solid line) and without control (dashed line); (b) phase portrait of the tip with control (solid line) and without control (dashed line); (c) error variatione1=x1−x˜1 and e2=x2−x˜2.
(a)(b)(c)The optimal linear feedback control by the SDRE is demonstrated to be effective in leading and maintaining the system in the desired orbit. Comparing the results obtained with the OLFC and SDRE controls, the results similarity is evident. The reason for that is because both controls are obtained considering an optimal control strategy.To analyze these controls in more detail, the next section presents a study of the parametric sensitivity of the controls.


### 3.3. Controlled System in the Presence of Parametric Errors

The control design is usually based on parameters of the mathematical model obtained from physical laws governed by the dynamic behavior of the system. Often, because of the limitations of the knowledge process, the used models do not accurately represent the real dynamics. Consequently, the control design cannot operate as intended when applied in a real process, because the parameters used in the control may contain parametric errors. To solve this problem, many researches have focused on incorporating the uncertainties associated with real structures into numerical simulation for reliable predictions [29]. The parametric uncertainties are associated with the discrepancies between the values of actual physical systems and the input parameters used for the analysis [29–31].To consider the effects of parameter uncertainties on the performance of the controller, the parameters used for the controls are considered as a random error of±20%; a similar strategy to that was used in [28], given by r^=r0.8+0.4ψτ, b^=b0.8+0.4ψτ, c^=c0.8+0.4ψτ, e^=e0.8+0.4ψτ, d^=d0.8+0.4ψτ, p^=p0.8+0.4ψτ, β^=β0.8+0.4ψτ, and Γ^=Γ0.8+0.4ψτ, where ψτ are normally distributed random functions. Hence, the analysis of the robustness of the controls is carried out in two ways: first, the random parameters are introduced only to the feedback controls; second, feedforward is introduced along with the feedback control with the uncertainties.Figure10 shows the deviation errors of the desired trajectory only for the feedback OLFC control considering the random parameters. The robustness of the control in keeping the system in the same orbit obtained with the control with uncertainties is observed, considering only the feedback control by OLFC (u), with the matrix of uncertainties given by(36)A=01−b^−r^Figure 10
Deviation of the desired trajectory with the proposed control with parametric errors only in feedback controlu by OLFC. (a) Deviation of the desired trajectory e1; (b) deviation of the desired trajectory e2.
(a)(b)Figures11(a) and 11(b) show the robustness of the SDRE control in keeping the system in the same desired orbit obtained accounting for the uncertainties only in the feedback control of the SDRE control (u) with the matrix A defined by(37)A=01−b^−c^e12+3e1x˜1+3x˜12−r^−p^/a^+e1+x˜13.Figure 11
Deviation of the desired trajectory with the proposed control with parametric errors only in feedback controlu by SDRE. (a) Deviation of the desired trajectory e1; (b) deviation of the desired trajectory e2.
(a)(b)In Figures12(a) and 12(b), the robustness of the control is in keeping the system in the desired orbit with the control considering the uncertainties in the feedback control by OLFC(u) and feedforward control u˜ with the matrix A denoted by (36) and(38)u˜=−x˜⋅2+r^x˜2+b^x˜1+c^x˜13−e^a^+x˜18+d^a^+x˜12+p^a^+x˜13x˜2−β^a^+x˜14−Γ^cosΩτ.Figure 12
Deviation of the desired trajectory with the proposed control with parametric errors in feedback control (u˜) by OLFC and feedforward control. (a) Deviation of the desired trajectory e1; (b) deviation of the desired trajectory e2.
(a)(b)In Figures13(a) and 13(b), we observed the robustness of the control in keeping the system in the desired orbit obtained along with the uncertainties in the feedback control by SDRE (u) and feedforward control u˜ with the matrix A(e) by (37) and(39)u˜=−x˜⋅2+r^x˜2+b^x˜1+c^x˜13−e^a^+x˜18+d^a^+x˜12+p^a^+x˜13x˜2−β^a^+x˜14−Γ^cosΩτ.Figure 13
Deviation of the desired trajectory with the proposed control with parametric errors in feedback controlu by SDRE and feedforward control u˜. (a) Deviation of the desired trajectory e1; (b) deviation of the desired trajectory e2.
(a)(b)Figure 14
Phase planes of the system in fractional-order derivative equations: (a)q3 = 1.0 and (b) q3 = 0.5031.
(a)(b)As can be seen in Figures10 and 11, the feedback control for both the OLFC control and the SDRE is robust to parametric variation. However, when the feedforward control signal is introduced (Figures 12 and 13), the feedforward control appears sensitive to parametric variation for both controls. Similar results are also observed in [28].


## 3.1. Time-Delayed Feedback Control (TDFC)

As originally suggested by the author in [25], a continuous control input to stabilize a chaotic oscillation is given by the difference between the current output and the past one as follows [12, 25–27]:(11)U=κgx1τ−T,x2τ−T−gx1τ,x2τ,where T is the time delay and κ is the feedback gain. The term gx1τ,x2τimplies a scalar output signal measured at the current time T and the past time τ−T, respectively. Since the control input equation (11) only depends on the output signal, the time delay T is adjusted to the period of a targeted stable periodic orbit that is intended to be stabilized in a chaotic attractor. The control input, therefore, converges to null after the controlled system is stabilized to the targeted orbit.Assuming that the velocity of oscillation is measured as an output of the nonlinear system (9), the control signal is given by(12)U=κx1τ−T,x2τ−T.The NEMS system with the control signal of (12) is expressed in the following way:(13)x˙1=x2,x˙2=−rx2−bx1−cx13+ea+x18−da+x12−pa+x13x2+βa+x14+ΓcosΩτ+U.The time delayT and feedback gain κ are important control parameters that substantially affect the control performance. The time delay T is adjusted to T=2π/Ω=2π/1=2π to stabilize an orbit with the same frequency as the external force oscillating the system. Figure 6 shows a bifurcation diagram of the system related to the gain κ.Figure 6
Bifurcation diagram for control gainκ.For anyκ≥0.06, the controlled system has a periodic behavior, with unitary period. Then, the gain is defined as κ=1.7 which is a value far from the chaotic region and that guarantees that the system will remain in the periodic orbit with unitary period. Figure 7 shows the comparison of the uncontrolled system (9) to the controlled one (13).Figure 7
Time histories of (a) displacement and (b) velocity; (c) phase portrait of the system with controlled trajectory (solid line) and the trajectory without control (dashed line); (d) control signalU.
(a)(b)(c)(d)The obtained periodic orbits with the time delay control of Figure7 can be numerically calculated through a Fourier series as(14)x1=−0.61128654−0.26404525sinτ,x2=−0.26404525cosτ.The time-delayed control was efficient in leading the system to a periodic orbit. As can be seen in the results, the control only used the control signal when necessary to take the system to a periodic orbit of the system, as observed in Figure7(d), where the control signal (U) tends to zero as the system stabilizes in the periodic orbit (Figures 7(a) and 7(b)). These results can also be observed in [12, 14].However, for the cases where it is desired to take the system to a previously defined orbit, the control may not be the most suitable. For these mentioned cases, the optimal control is considered in the next section.


## 3.2. Optimal Control

In the following, the solution of these problems using the optimal linear control technique for nonlinear systems developed by the authors in [28] is presented. The introduction of the control signal in (9) for the optimal control is considered equally in (13).However, the vector controlU for this optimal control consists of two parts: U=u˜+u, where u is the linear feedback control and u˜ is the feedforward control, and the control that maintains the system in the desired trajectory is given by(15)u˜=−x˜⋅2+rx˜2+bx˜1+cx˜13−ea+x˜18+da+x˜12+pa+x˜13x˜2−βa+x˜14−ΓcosΩτ.Substituting (15) into (13) and defining the deviation of the desired trajectory, (16)e=x1−x˜1x2−x˜2,where x˜1τand x˜2τ are the trajectories defined by (14).The system can be represented in the following form:(17)e˙1=e2,e˙2=−re2−be1−ce1+x˜13+ea+e1+x˜18−da+e1+x˜12−pa+e1+x˜13x˜2+βa+e1+x˜14+cx˜13−ea+x˜18+da+x˜12+pa+x˜13x˜2−βa+x˜14+u.


### 3.2.1. Feedback Control by OLFC

The optimal linear feedback control is applied as it was introduced by the authors in [28]. The system of (17) represented in deviations form can be written in matrix form as(18)e˙=Ae+Ge,x˜+Bu.According to [12], if there are matrices Q and R, positive and definite, with Q symmetric matrix, such that(19)Q˜=Q−GTe,x˜P−PGe,x˜,being positive and definite, in which the matrix G is limited, then the control u is optimal and it transfers the nonlinear system of (18) from any initial state to the final state e∞=0.Minimizing the cost functional to(20)J=∫0∞eTQ˜e+uTRudτ.The controlu can be found solving the following equation:(21)u=−R−1BTPe.HavingP as a symmetric matrix, the algebraic Riccati equation is developed, denoted by(22)PA+ATP−PBR−1BTP+Q=0.The matricesA and B may be represented by(23)A=01−b−r,B=01,and by definition(24)Q=10000100,R=1.Solving (22), the values of the cells of the matrix P are obtained, given by(25)P=109.49999.9019.90110.8503.Substituting the matricesR, B, and P into (21), the OLFC control signal, in deviation form, is given by(26)u=−9.9501e1−10.85504e2=−9.9501x1−x˜1−10.85504x2−x˜2.In Figure8, we observed the controlled system of (13) with control (U=u˜+u) in the designed orbit ((14).Figure 8
(a) Displacement of the tip with control (solid line) and without control (dashed line); (b) phase portrait of the tip with control (solid line) and without control (dashed line); (c) error variatione1=x1−x˜1 and e2=x2−x˜2.
(a)(b)(c)According to Figure8(c), the optimal linear feedback control was demonstrated to be effective in leading and maintaining the system in the desired orbit. Comparing the results obtained with the TDFC and OLFC controls, the OLFC control leads the system to the desired orbit much faster than the TDFC control. This is one of the advantages of using state feedback control (such as OLFC and SDRE). However, the determination of a desired orbit is needed, and, to keep the system in the desired orbit, the control signal has to be kept always acting in the system [12], which can be a disadvantage depending on the application of the controller. Nevertheless, a more suitable periodic orbit could be chosen depending on the purpose of the applied control.


### 3.2.2. Feedback Control by SDRE

For the SDRE control, (17) may be represented in matrix form as(27)e˙=Aee+Ge,x˜+Bu,where the matrix of states A(e) also depends on the deviations e.Considering that, the matricesA and B may be represented by(28)A=01−b−ce12+3e1x˜1+3x˜12−r−p/a+e1+x˜13,B=01.In addition, by definition(29)Q=10000100,R=1.The quadratic performance measured for the feedback control problem is given by(30)J=12∫t0∞eTQe+uTRudτ,where Q and R are positive definite matrices. Assuming full state feedback, the control law is given by(31)uf=−R−1BTPee=−kee.HavingPe as a symmetric matrix, the algebraic Riccati equation is denoted by(32)PeAe+AeTPe−PeQR−1BTPe+Q=0,keis obtained for each iteration by solving Riccati equation (32), and, considering the asymptotically stable periodic orbit in (14) as desired orbits (xj=xj∗), the optimal (or suboptimal) control problem (u) can be formulated as follows: determine the control signal uf ((31)) that transfers (9) from the initial state(33)x10=x10,x20=x20,into the final state(34)x1∞=x˜1,x2∞=x˜2,and minimizes the functional (30). Such minimization shown in (30) implies the minimization of the system deviation ((27)) of the desired state (e=x−x˜) and of the state of the control applied to the feedback control (u).Another important factor to consider is that the matrixA(e) cannot violate the controllability of the system. The system of (13) is controllable if the rank of the matrix M is 2, according to the rule(35)M=B2×1A2×2eB2×1.Then, to obtain a suboptimal solution for the dynamic control problem, the SDRE technique has the following procedure:(1)
Define the state-space model with the state-dependent coefficient as in (27).(2)
Definex(0) = x0, so that the rank of M is n and choose the coefficients of weight matrices Q and R.(3)
Solve Riccati equation (32) for the state xτ.(4)
Calculate the input signal from (31).(5)
Integrate system (27) and update the state of the system xτ with these results.(6)
Calculate the rank of (35); if rank = 2, go to step 3. However, if rank <2, the matrix A(e) is not controllable; therefore, the last matrix controllable A(e) that has been obtained should be used, and thus go to step 3.Figure9 shows system (13) without control and with control (U=u˜+u) in designed orbit equation (14).Figure 9
(a) Displacement of the tip with control (solid line) and without control (dashed line); (b) phase portrait of the tip with control (solid line) and without control (dashed line); (c) error variatione1=x1−x˜1 and e2=x2−x˜2.
(a)(b)(c)The optimal linear feedback control by the SDRE is demonstrated to be effective in leading and maintaining the system in the desired orbit. Comparing the results obtained with the OLFC and SDRE controls, the results similarity is evident. The reason for that is because both controls are obtained considering an optimal control strategy.To analyze these controls in more detail, the next section presents a study of the parametric sensitivity of the controls.


## 3.2.1. Feedback Control by OLFC

The optimal linear feedback control is applied as it was introduced by the authors in [28]. The system of (17) represented in deviations form can be written in matrix form as(18)e˙=Ae+Ge,x˜+Bu.According to [12], if there are matrices Q and R, positive and definite, with Q symmetric matrix, such that(19)Q˜=Q−GTe,x˜P−PGe,x˜,being positive and definite, in which the matrix G is limited, then the control u is optimal and it transfers the nonlinear system of (18) from any initial state to the final state e∞=0.Minimizing the cost functional to(20)J=∫0∞eTQ˜e+uTRudτ.The controlu can be found solving the following equation:(21)u=−R−1BTPe.HavingP as a symmetric matrix, the algebraic Riccati equation is developed, denoted by(22)PA+ATP−PBR−1BTP+Q=0.The matricesA and B may be represented by(23)A=01−b−r,B=01,and by definition(24)Q=10000100,R=1.Solving (22), the values of the cells of the matrix P are obtained, given by(25)P=109.49999.9019.90110.8503.Substituting the matricesR, B, and P into (21), the OLFC control signal, in deviation form, is given by(26)u=−9.9501e1−10.85504e2=−9.9501x1−x˜1−10.85504x2−x˜2.In Figure8, we observed the controlled system of (13) with control (U=u˜+u) in the designed orbit ((14).Figure 8
(a) Displacement of the tip with control (solid line) and without control (dashed line); (b) phase portrait of the tip with control (solid line) and without control (dashed line); (c) error variatione1=x1−x˜1 and e2=x2−x˜2.
(a)(b)(c)According to Figure8(c), the optimal linear feedback control was demonstrated to be effective in leading and maintaining the system in the desired orbit. Comparing the results obtained with the TDFC and OLFC controls, the OLFC control leads the system to the desired orbit much faster than the TDFC control. This is one of the advantages of using state feedback control (such as OLFC and SDRE). However, the determination of a desired orbit is needed, and, to keep the system in the desired orbit, the control signal has to be kept always acting in the system [12], which can be a disadvantage depending on the application of the controller. Nevertheless, a more suitable periodic orbit could be chosen depending on the purpose of the applied control.


## 3.2.2. Feedback Control by SDRE

For the SDRE control, (17) may be represented in matrix form as(27)e˙=Aee+Ge,x˜+Bu,where the matrix of states A(e) also depends on the deviations e.Considering that, the matricesA and B may be represented by(28)A=01−b−ce12+3e1x˜1+3x˜12−r−p/a+e1+x˜13,B=01.In addition, by definition(29)Q=10000100,R=1.The quadratic performance measured for the feedback control problem is given by(30)J=12∫t0∞eTQe+uTRudτ,where Q and R are positive definite matrices. Assuming full state feedback, the control law is given by(31)uf=−R−1BTPee=−kee.HavingPe as a symmetric matrix, the algebraic Riccati equation is denoted by(32)PeAe+AeTPe−PeQR−1BTPe+Q=0,keis obtained for each iteration by solving Riccati equation (32), and, considering the asymptotically stable periodic orbit in (14) as desired orbits (xj=xj∗), the optimal (or suboptimal) control problem (u) can be formulated as follows: determine the control signal uf ((31)) that transfers (9) from the initial state(33)x10=x10,x20=x20,into the final state(34)x1∞=x˜1,x2∞=x˜2,and minimizes the functional (30). Such minimization shown in (30) implies the minimization of the system deviation ((27)) of the desired state (e=x−x˜) and of the state of the control applied to the feedback control (u).Another important factor to consider is that the matrixA(e) cannot violate the controllability of the system. The system of (13) is controllable if the rank of the matrix M is 2, according to the rule(35)M=B2×1A2×2eB2×1.Then, to obtain a suboptimal solution for the dynamic control problem, the SDRE technique has the following procedure:(1)
Define the state-space model with the state-dependent coefficient as in (27).(2)
Definex(0) = x0, so that the rank of M is n and choose the coefficients of weight matrices Q and R.(3)
Solve Riccati equation (32) for the state xτ.(4)
Calculate the input signal from (31).(5)
Integrate system (27) and update the state of the system xτ with these results.(6)
Calculate the rank of (35); if rank = 2, go to step 3. However, if rank <2, the matrix A(e) is not controllable; therefore, the last matrix controllable A(e) that has been obtained should be used, and thus go to step 3.Figure9 shows system (13) without control and with control (U=u˜+u) in designed orbit equation (14).Figure 9
(a) Displacement of the tip with control (solid line) and without control (dashed line); (b) phase portrait of the tip with control (solid line) and without control (dashed line); (c) error variatione1=x1−x˜1 and e2=x2−x˜2.
(a)(b)(c)The optimal linear feedback control by the SDRE is demonstrated to be effective in leading and maintaining the system in the desired orbit. Comparing the results obtained with the OLFC and SDRE controls, the results similarity is evident. The reason for that is because both controls are obtained considering an optimal control strategy.To analyze these controls in more detail, the next section presents a study of the parametric sensitivity of the controls.


## 3.3. Controlled System in the Presence of Parametric Errors

The control design is usually based on parameters of the mathematical model obtained from physical laws governed by the dynamic behavior of the system. Often, because of the limitations of the knowledge process, the used models do not accurately represent the real dynamics. Consequently, the control design cannot operate as intended when applied in a real process, because the parameters used in the control may contain parametric errors. To solve this problem, many researches have focused on incorporating the uncertainties associated with real structures into numerical simulation for reliable predictions [29]. The parametric uncertainties are associated with the discrepancies between the values of actual physical systems and the input parameters used for the analysis [29–31].To consider the effects of parameter uncertainties on the performance of the controller, the parameters used for the controls are considered as a random error of±20%; a similar strategy to that was used in [28], given by r^=r0.8+0.4ψτ, b^=b0.8+0.4ψτ, c^=c0.8+0.4ψτ, e^=e0.8+0.4ψτ, d^=d0.8+0.4ψτ, p^=p0.8+0.4ψτ, β^=β0.8+0.4ψτ, and Γ^=Γ0.8+0.4ψτ, where ψτ are normally distributed random functions. Hence, the analysis of the robustness of the controls is carried out in two ways: first, the random parameters are introduced only to the feedback controls; second, feedforward is introduced along with the feedback control with the uncertainties.Figure10 shows the deviation errors of the desired trajectory only for the feedback OLFC control considering the random parameters. The robustness of the control in keeping the system in the same orbit obtained with the control with uncertainties is observed, considering only the feedback control by OLFC (u), with the matrix of uncertainties given by(36)A=01−b^−r^Figure 10
Deviation of the desired trajectory with the proposed control with parametric errors only in feedback controlu by OLFC. (a) Deviation of the desired trajectory e1; (b) deviation of the desired trajectory e2.
(a)(b)Figures11(a) and 11(b) show the robustness of the SDRE control in keeping the system in the same desired orbit obtained accounting for the uncertainties only in the feedback control of the SDRE control (u) with the matrix A defined by(37)A=01−b^−c^e12+3e1x˜1+3x˜12−r^−p^/a^+e1+x˜13.Figure 11
Deviation of the desired trajectory with the proposed control with parametric errors only in feedback controlu by SDRE. (a) Deviation of the desired trajectory e1; (b) deviation of the desired trajectory e2.
(a)(b)In Figures12(a) and 12(b), the robustness of the control is in keeping the system in the desired orbit with the control considering the uncertainties in the feedback control by OLFC(u) and feedforward control u˜ with the matrix A denoted by (36) and(38)u˜=−x˜⋅2+r^x˜2+b^x˜1+c^x˜13−e^a^+x˜18+d^a^+x˜12+p^a^+x˜13x˜2−β^a^+x˜14−Γ^cosΩτ.Figure 12
Deviation of the desired trajectory with the proposed control with parametric errors in feedback control (u˜) by OLFC and feedforward control. (a) Deviation of the desired trajectory e1; (b) deviation of the desired trajectory e2.
(a)(b)In Figures13(a) and 13(b), we observed the robustness of the control in keeping the system in the desired orbit obtained along with the uncertainties in the feedback control by SDRE (u) and feedforward control u˜ with the matrix A(e) by (37) and(39)u˜=−x˜⋅2+r^x˜2+b^x˜1+c^x˜13−e^a^+x˜18+d^a^+x˜12+p^a^+x˜13x˜2−β^a^+x˜14−Γ^cosΩτ.Figure 13
Deviation of the desired trajectory with the proposed control with parametric errors in feedback controlu by SDRE and feedforward control u˜. (a) Deviation of the desired trajectory e1; (b) deviation of the desired trajectory e2.
(a)(b)Figure 14
Phase planes of the system in fractional-order derivative equations: (a)q3 = 1.0 and (b) q3 = 0.5031.
(a)(b)As can be seen in Figures10 and 11, the feedback control for both the OLFC control and the SDRE is robust to parametric variation. However, when the feedforward control signal is introduced (Figures 12 and 13), the feedforward control appears sensitive to parametric variation for both controls. Similar results are also observed in [28].


## 4. AFM-TM Modeling with Fractional-Order Differential Equation

Due to the small scale of the AFM operation, a viscoelastic behavior of the system is observed. For the viscoelastic behavior, it is possible to use the fractional-order derivatives to analyze the dynamic behavior of the system [31]. Due to the essential differences between ordinary differential equations (ODE) and fractional-order differential equations (FODEs), most of the characteristics or conclusions of the ODE systems cannot be directly extended to the case of the FODE systems. Differential equations may involve Riemann-Liouville differential operators of fractional order q > 0, which generally takes the form [32–34](40)Dqxt=1Γη−q∫τ0τxηuτ−uq−η+1du,where η is the first integer not less than q. It is easily proved that the definition is the usual derivatives definition when q = 1. The case for 0 < q < 1 seems to be particularly important [10]. For simplicity and without loss of generality, in the following, it is assumed that t0 = 0, 0 < q < 1. The FODEs are numerically integrated using the algorithm proposed by Petras [33] with integration step of h = 0.0001.Hence, the technique of fractional calculus is included to analyze the behavior of system (9) with the squeeze-film damping −px2/a+x13 as a fractional-order term, yielding(41)dq1x1dτq1=x2,dq2x2dτq2=x2=−rx2−bx1−cx13+ea+x18−da+x12−pa+x13dq3x3dτq3+βa+x14+ΓcosΩτ,dq3x3dτq3=x2,where 0<q1,q2,q3≤1and the fractional order is denoted by q=q1,q2,q3.Figures15(a) and 15(b) show the phase plane of (41) for the derivative orders as q3 = 1.0 and q3 = 0.5031, respectively. In addition, Figure 15 shows the time history of displacement for both q3 = 1.0 and q3 = 0.5031. Those results show an evidence of irregular motion, which is an evidence of chaotic behavior.Figure 15
Time series of fractional order of the system withq3 = 1.0 and q3 = 0.5031.To calculate the Lyapunov exponent for system (9), a second-order ODE is considered as in the method proposed by [24], which uses the eigenvalues of the Jacobian matrix of classical ODE system (9). However, the fractional-order derivative does not allow the same method due to the fractional derivative. Therefore, 0-1 test is considered to analyze whether the system is chaotic or not. The 0-1 test has proven to be a very useful and practical numerical tool to identify the behavior of dynamical systems and can be found in similar problems that were used with high success [14, 21].


### 4.1. The 0-1 Test

The 0-1 test, proposed by the authors in [35–37], is directly applied to a time series data based on the statistical properties of a single coordinate, here in the variable x1 of (41). Basically, the 0-1 test consists of estimating a single parameter Kc. The test considers a system variable xj, where two new coordinates (p,q) are defined as follows:(42)pn,c¯=∑j=0nxjcosjc¯,(43)qn,c¯=∑j=0nxjsinjc¯,where c¯∈0,π is a constant. The mean square displacement of the new variables pn,c¯ and qn,c¯ is described by(44)Mn,c=limn⟶∞1N∑j=1Npj+n,c¯−pj,c¯2+qj+n,c¯−qj,c¯2,where n=1,2,…,N and, therefore, we obtain the parameter Kc in the limit of a long time Kc=covY,Mc¯/varYvarMc¯ through(45)Kc=covY,Mc¯varYvarMc¯,where vectors Mc¯=M1,c¯,M2,c¯,K,Mnmax,c¯ and Y=1,2,…,nmax.Given any two vectorsx and y, the covariance covx,y and variance varx, of nmax elements, are usually defined as [31](46)covx,y=1nmax∑n=1nmaxxn−x¯yn−y¯,varx=covx,x,where x¯ and y¯ are the average of xn and yn, respectively. The value of the parameter K is obtained by (45). If the Kcvalue is close to 0, the system is periodic; on the other hand, if the Kc value is close to 1, the system is chaotic.Figure16 shows the values of the parameter Kc of the 0-1 test for a scan of the fractional derivative q3∈0.5:1 of the AFM-TM system to verify the presence of the chaotic behavior of the system. As it is observed, the system has a chaotic behavior for all values of q3 mainly to q3=0.5130 which was previously analyzed and whose Kc is very close to 1.Figure 16
Representation of Kc (0-1 test) vs.0.5≤q3≤1 (fractional order).As forq3=0.5130the system presents chaotic behavior; this derivative value is adopted to the control analysis of the system in fractional order.


## 4.1. The 0-1 Test

The 0-1 test, proposed by the authors in [35–37], is directly applied to a time series data based on the statistical properties of a single coordinate, here in the variable x1 of (41). Basically, the 0-1 test consists of estimating a single parameter Kc. The test considers a system variable xj, where two new coordinates (p,q) are defined as follows:(42)pn,c¯=∑j=0nxjcosjc¯,(43)qn,c¯=∑j=0nxjsinjc¯,where c¯∈0,π is a constant. The mean square displacement of the new variables pn,c¯ and qn,c¯ is described by(44)Mn,c=limn⟶∞1N∑j=1Npj+n,c¯−pj,c¯2+qj+n,c¯−qj,c¯2,where n=1,2,…,N and, therefore, we obtain the parameter Kc in the limit of a long time Kc=covY,Mc¯/varYvarMc¯ through(45)Kc=covY,Mc¯varYvarMc¯,where vectors Mc¯=M1,c¯,M2,c¯,K,Mnmax,c¯ and Y=1,2,…,nmax.Given any two vectorsx and y, the covariance covx,y and variance varx, of nmax elements, are usually defined as [31](46)covx,y=1nmax∑n=1nmaxxn−x¯yn−y¯,varx=covx,x,where x¯ and y¯ are the average of xn and yn, respectively. The value of the parameter K is obtained by (45). If the Kcvalue is close to 0, the system is periodic; on the other hand, if the Kc value is close to 1, the system is chaotic.Figure16 shows the values of the parameter Kc of the 0-1 test for a scan of the fractional derivative q3∈0.5:1 of the AFM-TM system to verify the presence of the chaotic behavior of the system. As it is observed, the system has a chaotic behavior for all values of q3 mainly to q3=0.5130 which was previously analyzed and whose Kc is very close to 1.Figure 16
Representation of Kc (0-1 test) vs.0.5≤q3≤1 (fractional order).As forq3=0.5130the system presents chaotic behavior; this derivative value is adopted to the control analysis of the system in fractional order.


## 5. Determination of the Optimal Control for Fractional-Order Case

In this section, the techniques of fractional calculus to analyze the behavior of system (41) due to the influence of the squeeze-film damping −px2/a+x13 are considered. Thus, a control signal is introduced into (41), according to(47)dq1x1dτq1=x2,dq2x2dτq2=x2=−rx2−bx1−cx13+ea+x18−da+x12−pa+x13dq3x3dτq3+βa+x14+ΓcosΩτ+U,dq3x3dτq3=x2.The vector controlU consists of two parts: U=u˜+u, where u˜ is the feedforward control and u is the linear feedback control. The defined periodic orbit is x˜1τx˜2τx˜3τT. If the function x˜1τx˜2τx˜3τTis the solution of (44), without the control U, then u = 0.In this way, the desired regime is obtained by(48)dq1x˜1dτq1=x˜2,dq1x˜2dτq2=−rx˜2−bx˜1−cx˜13+ea+x˜18−da+x˜12−pa+x˜13dq3x˜3dτq3+βa+x˜14+ΓcosΩτ+u˜,dq3x˜3dτq3=x˜2.The feedforward controlu˜ is given by(49)u˜=−dq1x˜2dτq2+rx˜2+bx˜1+cx˜13−ea+x˜18+da+x˜12+pa+x˜13dq3x˜3dτq3−βa+x˜14−ΓcosΩτ.Substituting (49) into (47) and defining the deviation of the desired trajectory(50)e=x1−x˜1x2−x˜2x3−x˜3.The system can be represented in the following form:(51)dq1e1dτq1=e2,dq2e2dτq2=−re2−be1−ce1+x˜13+ea+e1+x˜18−da+e1+x˜12−pa+e1+x˜13dq3e3+x˜3dτq3+βa+e1+x˜14+cx˜13−ea+x˜18+da+x˜12+pa+x˜13dq3e3dτq3−βa+x˜14+u,dq3e3dτq3=e2.Since the objective of this work is to controlx1 and x2, the variable x3 is considered only as a disturbance of the system, as similarly proposed in [14, 28].


### 5.1. Feedback by OLFC Control

The matricesA, B, Q, and R are given by(52)A=01−b^−r^,B=01,Q=10000100,R=1.Figure17 shows controlled system (47) with control signal U=u˜+u for the desired orbit of (14). It is observed that, even with the introduction of the fractional-order derivative, the OLFC control is efficient to control the system.Figure 17
OLFC control: (a) phase portrait with control (solid line) and without control (dashed line); (b) trajectory of the control (dashed line) in comparison with the desired trajectory (solid line); (c) error variatione1=x1−x˜1 and e2=x2−x˜2, by considering q1=q2=1 and q3=0.5031.
(a)(b)(c)


### 5.2. Feedback by SDRE Control

The SDRE control considers the matricesA, B, Q, and R as(53)A=01−b−ce12+3e1x˜1+3x˜12−r−p/a+e1+x˜13,B=01,Q=10000100,R=1.Figure18 shows controlled system (47) with control signal U=u˜+u for the desired orbit of (14). Similarly, the SDRE control is efficient to control the system in fractional-order derivative.Figure 18
SDRE control: (a) phase portrait with control (solid line) and without control (dashed line); (b) trajectory of the control (dashed line) in comparison with the desired trajectory (solid line); (c) error variatione1=x1−x˜1 and e2=x2−x˜2, by considering q1=q2=1 and q3=0.5031.
(a)(b)(c)As observed in Figures17(c) and 18(c), the optimal linear feedback control obtained by OLFC and SDRE control was demonstrated to be effective in leading and maintaining the system with squeeze-film damping and fractional order to the desired orbit.Both OLFC and SDRE control techniques are very efficient in controlling the chaotic AFM system which presented chaotic behavior due to the presence of squeeze-film damping considered as a viscous damping ((13)) and also viscoelastic damping ((51)). The viscoelastic damping is a complex contribution due to the need of the introduction of the fractional-order derivative, in which the SDRE and OLFC controls can dominate.


## 5.1. Feedback by OLFC Control

The matricesA, B, Q, and R are given by(52)A=01−b^−r^,B=01,Q=10000100,R=1.Figure17 shows controlled system (47) with control signal U=u˜+u for the desired orbit of (14). It is observed that, even with the introduction of the fractional-order derivative, the OLFC control is efficient to control the system.Figure 17
OLFC control: (a) phase portrait with control (solid line) and without control (dashed line); (b) trajectory of the control (dashed line) in comparison with the desired trajectory (solid line); (c) error variatione1=x1−x˜1 and e2=x2−x˜2, by considering q1=q2=1 and q3=0.5031.
(a)(b)(c)


## 5.2. Feedback by SDRE Control

The SDRE control considers the matricesA, B, Q, and R as(53)A=01−b−ce12+3e1x˜1+3x˜12−r−p/a+e1+x˜13,B=01,Q=10000100,R=1.Figure18 shows controlled system (47) with control signal U=u˜+u for the desired orbit of (14). Similarly, the SDRE control is efficient to control the system in fractional-order derivative.Figure 18
SDRE control: (a) phase portrait with control (solid line) and without control (dashed line); (b) trajectory of the control (dashed line) in comparison with the desired trajectory (solid line); (c) error variatione1=x1−x˜1 and e2=x2−x˜2, by considering q1=q2=1 and q3=0.5031.
(a)(b)(c)As observed in Figures17(c) and 18(c), the optimal linear feedback control obtained by OLFC and SDRE control was demonstrated to be effective in leading and maintaining the system with squeeze-film damping and fractional order to the desired orbit.Both OLFC and SDRE control techniques are very efficient in controlling the chaotic AFM system which presented chaotic behavior due to the presence of squeeze-film damping considered as a viscous damping ((13)) and also viscoelastic damping ((51)). The viscoelastic damping is a complex contribution due to the need of the introduction of the fractional-order derivative, in which the SDRE and OLFC controls can dominate.


## 6. Conclusion

This work presented the dynamical analysis and control of an AFM-TM system with the addition of the Casimir force and VdW forces that induced the presence of a chaotic behavior. In addition, the investigation of the fractional-order derivative is sought out to allude the influence of the viscoelastic term in the AFM-TM setup. A set of parameters enunciated the chaotic behavior. Thus, the intensity of the Casimir force strongly influenced the behavior of the AFM-TM, differently of the damping parameterp that only has a few sets of parameters that strongly influence the chaotic behavior.In order to suppress the chaotic motion, the TDFC and the optimal control by OLFC and SDRE techniques were projected. The efficiency of the proposed techniques was demonstrated through numerical simulations. As could be seen in Figure7(c), the time-delayed control led the system to one of the periodic attractors of the system. Thus, the control signal was necessary until the system stabilizes in the periodic orbit (Figure 7(d)). Using Fourier series, it is possible to determine the mathematical model of the periodic orbit obtained with the use of TDFC control.As the time-delayed control is not designed to take the system to any previously defined orbit, an alternative is the application of the OLFC or SDRE control. As could be seen in Figures9 and 10, the OLFC and the SDRE control were efficient in taking the system to the same periodic orbit obtained by the TDFC.Thus, it is possible to conclude that the TDFC is an excellent option when the objective is to take the system to a periodic orbit with the lowest cost of control. The OLFC or SDRE control is a viable option to be considered in cases in which it is necessary to impose the desired orbit. In addition, the robustness of the OLFC and SDRE controls was analyzed. Both OLFC and SDRE techniques worked well and were robust due to parametric error analyses.As the OLFC and SDRE control proved to be robust to parametric errors, the application of the control in the fractional-order system was also considered. Numerical results showed that the OLFC and SDRE controls are also effective for control in fractional-order systems. For the OLFC control, the feedback control is linear and not state dependent, so its processing is faster than the state-dependent SDRE feedback control. However, as it does not update the states in each step, the OLFC control is more sensitive to parametric errors in the feedforward control than the feedforward control used in the SDRE control, whose results are also observed in [27]. Thus, the OLFC control is the most suitable for the cases where the processing time is a limiter for control application.



---

*Source: 4048307-2020-06-22.xml*