Limit Cycles in Nonlinear Systems with Fractional Order Plants

In recent years, there has been considerable interest in the study of feedback systems containing processes whose dynamics are best described by fractional order derivatives. Various situations have been cited for describing heat flow and aspects of bioengineering, where such models are believed to be superior. In many situations these feedback systems are not linear and information on their stability and the possibility of the existence of limit cycles is required. This paper presents new results for determining limit cycles using the approximate describing function method and an exact method when the nonlinearity is a relay characteristic.


Introduction
In recent years, there has been considerable interest in the study of feedback systems containing processes whose dynamics are best described by fractional order derivatives.A physical system OPEN ACCESS represented by a differential equation where the orders of derivatives can take any real number and not necessarily an integer number can be called a fractional order system.The idea of non-integer order (or fractional order) differentiation/integration emerged in 1695 [1].At the beginning, mathematicians studied it only as a theoretical subject due to its complexity and therefore other science disciplines could not use it effectively owing to the absence of exact solution methods of non-integer order differential equations [2].However, in recent years, facilitated by today's computational facilities, considerable attention has been given to fractional order systems, including fractional order control systems.
Interest in fractional calculus has made a significant impact on the field of control engineering and as a result of this there have been many studies dealing with the analysis and design of control systems.It has been seen that the frequency domain approaches of classical control, where the Laplace complex variable of a transfer function is replaced by jω, such as Bode, Nyquist and Nichols diagrams can be applied to a fractional order plant without any modification.However, for time domain analysis integer approximations of fractional powers of s, such as given in [3][4][5][6], which are not exact have to be used.On the design side, fractional order versions of classical controllers, such as fractional order PID, have been designed [7] and applied [8][9][10][11] and various results have been obtained for stability analysis [12,13].However, the field of fractional order control systems still needs further research on many important problems.One of these situations is when nonlinearity occurs in a control system.It is well known that all practical systems are nonlinear and linear models of many engineering systems are always approximate [14].It is apparent, since as mentioned above the frequency domain representation of a fractional order plant is exact, that frequency domain methods can play an important role in analysis and design.
The aim of this paper is to study the problem of the stability of an autonomous nonlinear feedback system with a fractional order plant.The stability of a nonlinear control system can often be demonstrated by showing that there is no possibility of the existence of sustained oscillations, known as limit cycles.The most powerful, although approximate, method for such investigations is the Describing Function (DF) method [15,16].The frequency at the intersection point, when one exists, of the Nyquist plot of the system and the complex plot of the negative inverse of the DF is used in limit cycle analysis.A review of the classical DF method is given and application to a control system with a fractional order plant is demonstrated.Tsypkin [17] presented a frequency domain method for the exact analysis of limit cycles in relay feedback systems and the problem is formulated in terms of A loci [14], which are frequency dependent loci and more general than the Tsypkin loci.A program is developed in MATLAB for computation of limit cycles by the Tsypkin method in a feedback loop with an on-off relay, relay with hysteresis or relay with dead zone and a fractional order plant.
The rest of the paper is organized as follows: A review of fractional order dynamics is given in Section 2. Section 3 studies stability and limit cycles in nonlinear systems with FOTF.The DF analysis for nonlinear fractional order control systems is given.Computation of the A locus of a fractional order plant and some examples are presented for limit cycles in relay systems.Section 4 gives some concluding remarks.

Fractional Order Dynamics
Many real systems are known to display fractional order dynamics.For example, it is known that the semi-infinite lossy (RC) transmission line demonstrates fractional behavior since the current into the line is equal to the half derivative of the applied voltage that is V(s) = (1/√)I(s) [18].Thus, the significance of fractional order representation is that fractional order differential equations are more adequate to describe some real world systems than those of integer order models [19,20].Many physical systems such as viscoelastic materials [21,22], electromechanical processes [23], long transmission lines [24], dielectric polarisations [25], colored noise [26], cardiac behavior [27], problems in bioengineering [28], and chaos [29] can be described using fractional order differential equations.Thus, fractional calculus has been an important tool to be used in engineering, chemistry, physical, mechanical and other sciences.Extensive results on fractional order systems and control can be found in the book by Monje et al. [30].
Fractional order calculus is a generalization of the ordinary differentiations by non-integer derivatives.Many mathematicians like Liouville and Riemann contributed to the field of fractional calculus.There are different definitions of fractional order operators such as Grünwald-Letnikov, Riemann-Lioville and Caputo [31].The Caputo definition can be stated as [32], where D α y(t) = d α y(t)/dt α indicates the Caputo derivative of y(t), α ∈ R + is the rational number, [α] is the integer part of  and L denotes the Laplace transform.
A fractional order control system with input r(t) and output y(t) can be described by a fractional differential equation of the form [33], or by a fractional order transfer function of the form, where   ,   ( = 0,1,2, … ,  and  = 0,1,2, … ,  ) are real parameters and   ,   are real positive numbers and  0 <  1 < ⋯ <   and  0 <  1 < ⋯ <   .Thus, a transfer function including fractional powered s terms can be called a fractional order transfer function, FOTF.For example, with the FOTF () = 1 ( Bode and Nyquist diagrams of this equation can then be obtained as shown in Figures 1 and 2. This example is given to show that the frequency response computation of FOTF can be obtained similar to integer order transfer functions.However, for time domain computation, there is not a general analytical method for determining the output of a control system with an FOTF.There have been many studies over the years [34][35][36], some of them are based on integer approximation models and others based on numerical approximation of the non-integer order operator.The methods developing integer order approximations can be used for time domain analysis of fractional order control systems similar to classical control approaches.Some of the well known methods for evaluating rational approximations are the Continued Fractional Expansion (CFE) method, Oustaloup's method, Carlson's method, Matsuda's method, Chareff's method, least square methods and others [18,[37][38][39][40][41].In this paper, the CFE method given in [35,42] is used for time domain simulation of nonlinear control systems with FOTF.This method finds approximations from where  =  − 1.Thus, first, second, third and fourth integer order approximations of   can be found from the following equations The fourth order approximation has been used for the simulations reported in this paper.The following examples show that the fourth order approximations given in Equation (9) give good results.
The Bode plots of first, second, third, fourth order approximations of () =  0.5 and its exact Bode plots are shown in Figure 3 where it can be seen that the fourth order Bode plots fit with the exact Bode plots for the frequencies approximately from 0.04 rad/s to 20 rad/s.The step responses for () = 1/( 1.5 + 1) using first, second, third and fourth order approximations for  0.5 are given in Figure 4.The third and fourth order approximations match almost exactly.One of the most well-known approximation techniques is Oustaloup's method [37].Oustaloup's method gives 1,3,5,7,… order integer approximations.A comparison between approximations given in Equation ( 9) and Oustaloup's method for () = 1/( 1.5 + 1) is given in Figure 5.The comparison is done between the seventh order Oustaloup's approximation and the fourth order approximation given in Equation ( 9).The seventh order Oustaloup's approximation for  0.5 is  0.5 ≅ 10 7 + 509.4 6 + 5487 5 + 14990 4 + 10790 3 + 2045 2 + 98.34 + 1  7 + 98.34 6 + 2045 5 + 10790 4 + 14990 3 + 5487 2 + 509.4 + 10 (10) Figure 5 shows that two step responses agree very closely.However, we prefer to use the fourth order approximation given in Equation ( 9) since Oustaloup's approximation is seventh order.Also in [35] it has been shown that the approximations given in Equations ( 6)-( 9) provide very good results.Step responses for () = 1/( 1.5 + 1) using first, second, third and fourth order approximations given in Equations ( 6)- (9).Step responses using the fourth order approximation given in Equation ( 9) and seventh order Oustaloup approximation.

Describing Function Method
The DF, (), of a nonlinear element is defined as the ratio of the fundamental output to the magnitude of an applied sinusoidal input [43].Considering a nonlinear element n(x) with input  =  and corresponding output (), then if n(x) has odd symmetry, the fundamental component in (), namely  1  +  1 , has The describing function () is then given by which will be real, that is,  1 = 0, if the nonlinearity is single valued.To investigate the possibility of a limit cycle in the system of Figure 6 the characteristic equation is then examined.Typically, this is done using a Nyquist diagram where the loci () and () = −1 () ⁄ are plotted, and any intersection of the loci gives the amplitude and frequency of a possible limit cycle.Fortunately, the frequency domain analysis of a FOCS (Fractional Order Control System) can be conducted in a similar way to that of an integer order one.Therefore, the frequency domain expression can be easily obtained by substituting  =  in the Laplace transform of the transfer function.Since the describing function method is a frequency domain approach, it can be applied to the FOCS to analyze some aspects of the effect of nonlinearity on its performance.When the nonlinearity in the negative unity feedback control system of Figure 6 is an ideal relay or relay with hysteresis as shown in Figures 7 and 8 nonlinearity then their DFs are respectively: The possibility of a limit cycle can then be investigated as shown in the following example.
Example 1: Consider the control system of Figure 6 with the following fractional order plant   ( ) presents the Nyquist plot of the plant () and the negative inverse of the describing function in Equation (15).From Figure 9, a limit cycle of frequency  = 1.632 rad/s is predicted for the fractional order system.Since () = 4ℎ  ⁄ for the ideal relay nonlinearity, for ℎ = 1, the approximate amplitude of the limit cycle can be calculated from the intersection point of () and () as, This amplitude, according to DF theory, is an approximation for the amplitude of the fundamental, although it is often assumed for convenience to be the amplitude of the distorted limit cycle.The period of the oscillation is given by; The simulation results for the nonlinear control system with the relay nonlinearity are given in Figure 10.The simulation results, for which the measured peak amplitude and period were 1.025 and 3.93 s respectively, agree well with the DF analysis method although both are approximate.As the limit cycle is near to sinusoidal, the DF result may be expected to be quite accurate.

Tsypkin's Method
The analysis of the relay feedback problem can be studied using the frequency domain approach of the Tsypkin method and the A Loci approach [14].For a feedback loop containing a relay and a transfer function (), limit cycles can be found using the A loci of () where   (, ) =   (, ) +   (, ) with and where () =   () +   ().  (, ) is a generalized summed frequency locus with its real and imaginary values at a particular frequency, , depending on weighted values, according to the choice of, , of the real and imaginary values of () at the frequencies  for  = 1,2, … , ∞.In particular, for  = 0 the real (imaginary) part of   (, ) depends only upon the real (imaginary) part of () at frequencies .It is possible to evaluate   (, ) and   (, ) by computationally summing to a finite number, M, of terms.For integer transfer functions, closed form solutions have been found for the infinite series but this remains an open problem for fractional order transfer functions (FOTFs).When the limit cycle is odd symmetrical only odd values of n are required in the series and the corresponding A locus is denoted by   .Some properties of the A loci are: The   locus for  = 0 is identical with the Tsypkin locus [17], Λ(), if lim →∞ () = 0, apart from a constant factor.The relationship is The   locus can thus be regarded as a generalized Tsypkin locus.The A locus satisfies the superposition property, that is if the linear plant () =  1 () +  2 (), then If  1 () = () − , then The loci are periodic in , with period 2, that is For   the periodicity is odd, that is Since any integer plant transfer function can be written in terms of a summation of transfer functions having a real or complex pair of poles, use of Equation (24) allows A loci for most integer transfer functions to be obtained in terms of the A loci of a few basic transfer functions.
For a relay with no dead zone, the simplest assumption for the periodic output is a square wave with 1:1 mark space ratio.By taking the Fourier series for the square wave and calculating the corresponding series for the output of  it can be shown that the output and the derivative of the output of the linear transfer function  can be expressed in terms of the A loci [14].This also applies to the situation where the relay has dead zone, where unlike the situation for the relay without dead zone, finite values of  are required.Taking the input to the relay, () = −(), where c(t) is the output of G(s), and assuming the relay to have a hysteresis of ±∆, that is it switches from −h to +h for an input of ∆ at time zero and from +h to −h for an input of −∆, the switching condition (0) = ∆ and (0) > 0 (28) yields the result    (0, ) = − Δ 4ℎ and  ⁄    (0, ) < 0 provided lim →∞ () = 0.When the latter is not the case, discontinuities in the relay input or its derivative may occur at the switching instants and known corrections need to be applied to the equations.For a relay with dead zone the solution is obtained from (0, ) −    (Δ, ) < 0 and (0, ) −    (−Δ, ) < 0 Here, there are two equations as there are two unknowns the pulse width, Δ, and the limit cycle frequency, .Note the expressions involving the real parts need not normally be checked but it is easily done if the solutions are obtained graphically by plotting for the relay with no dead zone    (0, ) and for the relay with dead zone    (0, ) −    (Δ, ) and    (0, ) −    (−Δ, ).
Consider the transfer function () = 1 ( + 1)( + 2) ⁄ in a feedback loop having a relay with hysteresis.For this transfer function thus for an odd symmetrical limit cycle The limit cycle solution frequency is given where the above locus meets the line − Δ 4ℎ ⁄ parallel to the negative real axis on a Nyquist plot.The DF solution is given where () meets this line, so when a program is written where  can be input then the DF solution can be obtained for  = 1 and convergence can be seen by taking larger and larger values of .Thus, plotting    (0, ) and () provides a good graphical approach for obtaining the limit cycle solution.Another option, to obtain the limit cycle frequency, is to solve for .Strictly speaking, one should also check that the real part of the locus is negative but this is obviously the case for this simple example.Cases where it is not typically involve "unusual" transfer functions, and the DF solution may then also be a problem.For the relay with dead zone, two sets of loci must be plotted for the graphical approach because the two loci must be plotted for different values of  =  say, and the solution is obtained when the given lines are intersected with the same  and  on the loci.Figure 11 shows A Loci of Equation ( 34) for different values of n.Table 1 shows the solution of Equation ( 36) for the relay with hysteresis for different values of n.Convergence is fast because the transfer function is a good low pass filter.Example 2: Consider the nonlinear system of Figure 6 with the integer order transfer function  2 () = 2 −  ( + 2) 2 (37) and an ideal relay.This transfer function is not a good filter so there is some distance between the A locus and Nyquist plots as seen in Figure 12.The frequency of the limit cycle using the A locus plot, as seen from Figure 12 is  = 2.8985 rad/s for n=101 and  = 3.4641 rad/s from the Nyquist plot.The limit cycle is shown in Figure 13, which can be seen to be quite distorted, hence the quite large difference in the two frequencies.Solving the equation corresponding to Equation (36) for this example, for different values of n, gives Table 2.This shows, as expected, that larger values of n must be used to get the limit cycle frequency accurately.A Loci for different values of n are shown in Figure 14.(37) where () is an ideal relay.
Example 3: In this example, the Typskin loci for two plants, one is fractional order, with time delay are studied.These are: and The A Locus and Nyquist plots for  3 () and  3 () are shown in Figures 15 and 16, respectively.
From the figures, the exact limit cycle frequencies are 1.27 rads/s and 0.86 rads/s, respectively, for the integer order and fractional order plants.Example 4: Consider Figure 6 with the fractional order plant  4 () = 2  3.6 + 3 2.4 + 3 1.2 + 1 (40) Assume that the nonlinear system includes a relay with no dead zone.The A Locus, and the () loci for the ideal relay and relay with hysteresis, having ℎ Δ = 2 ⁄ , are shown in Figure 17. Figure 18 shows the time responses for the ideal relay and relay with hysteresis.The limit cycle frequencies are  = 1.1321 rad/s and  = 0.9681 rad/s, respectively.
It can be shown using the approximate DF theory that a necessary but not sufficient criterion for stability of a limit cycle is that the intersection between () and () should be as shown in Figure 15.That is, if one moves along () in the direction of increasing frequency then at the point of intersection increasing amplitude on the () locus will lie to the left.
in a feedback loop having a relay with dead zone as shown in Figure 19.The DF for a relay with dead zone only, that is for  > , Δ = 0 is and its input and output waveforms are shown in Figure 20.The real and imaginary parts of the denominator of  5 () are −0.309ω 1.2 -1.902ω 2.2 + 0.309ω 3.2  and 0.9510ω 1.2 -0.6180ω 2.2 -0.9510ω 3.2 , respectively.Thus, Equation (14), which is 1 + N(a)G(jω) = 0 gives Figure 22 shows a plot of () and () from which it is seen that the latter travels along the negative real axis from minus infinity as a increases to a maximum of − 2ℎ ⁄ and then returns.Thus, there are two intersections with () and according to the above criterion only the larger amplitude one corresponds to a stable limit cycle.

Analysis of Limit Cycle Existence According to Compensator Gain
The above transfer function for K = 1 crosses the negative real axis at a frequency of 0.7265 rad/s and its magnitude is given by: 1 () 1.2 An intersection with () will thus only occur if Taking  = 1, ℎ =  for the relay gives ()  = 2, so that for stability K < 0.521.Figure 23 shows simulation results for K = 0.49 and 0.55.The former is seen to be stable and a limit cycle exists for the latter.In the simulation, a limit cycle developed for K just greater than 0.49.Thus, although both results are approximate, there is good agreement.In this example, the effect of the change in order of a fractional integrator on the feedback loop stability is investigated.
Example 6: Consider the control system of Figure 19 with the fractional order plant   () = 2   ( + 1) 2  (47 where  is the fractional order of the integrator 1   ⁄ .The relay is again taken to have  = 1, ℎ = .The () loci for  = 1.2 and  = 0.779, which passes through −0.5, are shown in Figure 24.Example 7: Here, example 5 shown in Figure 19 is again considered to illustrate the Tsypkin method.The relay parameters were again taken as  = 1, ℎ =  to allow comparison of the results with those obtained using the DF method and initially K was taken equal to 1.The required A loci plots    (0, ) −    (Δ, ) and    (0, ) −    (−Δ, ) were plotted for a selection of values of ∆ as shown in Figure 26, and the values of  and  were recorded where they met the lines − 2ℎ ⁄ and  2ℎ ⁄ , respectively.Table 4 shows recorded values of  and  for K = 1.    4.
The procedure was repeated for decreasing values of K and the results are given in Table 5 which shows the frequency and pulse width of the stable and unstable limit cycles.There are two frequency values for each K and larger one corresponds to a stable limit cycle.The value of K for which no solution was possible, and thus the system was stable, was 0.49.The results obtained from simulation, DF and A locus methods for comparison are given in Table 6 where it can be seen that the system is stable for  ≤ 0.52 according to DF method and it is stable for  ≤ 0.49 according to simulation and A locus methods.Time domain simulation using Figure 19 is shown in Figure 28 which is the stable limit cycle obtained from a simulation with K = 1.Limit cycle frequency and pulse width are computed as  = 0.7173 rad/s and  = 3.65, respectively.

Conclusions
In this paper, after a brief review of the DF method the computation of the A locus for a fractional order plant has been studied.The A locus method is important as it allows calculation of the exact limit cycle frequency in a relay control system.For the fractional order plant, the A locus has to be computed by summing terms of a series and taking one term only is the equivalent of the DF method.Several examples have been given showing applications of the DF and A loci methods to the computation of limit cycles in fractional order plants and the results compared.It is known that the calculations to find the limit cycles in feedback systems containing relay type characteristics can be done using either time domain or frequency domain methods.Therefore, the problem of the calculation of limit cycles in a system with a relay with dead zone and a fractional order transfer function using time domain approach can be considered as future work and comparing the exact frequency domain approach with the approximate time domain solution can provide useful results.

Figure 3 .
Figure 3. Exact Bode plots of  0.5 and first, second, third, fourth order approximations.

Figure 5 .
Figure 5.Step responses using the fourth order approximation given in Equation (9) and seventh order Oustaloup approximation.

Figure 7 .
Figure 7.The diagram of the ideal relay.

Figure 8 .
Figure 8.The diagram of the relay with hysteresis.

Figure 10 .
Figure10.Time responses of closed loop system in Figure6with  1 () and ideal relay.

Figure 19 .
Figure 19.A nonlinear feedback system having relay with dead zone.

Figure 23 .
Figure 23.Time responses of closed loop system in Figure 19 with Equation (43) and relay with dead zone for K = 0.49 and 0.55

Figure 25 .
Figure 25.Nyquist plots for different values of  for Equation (47) where () is a relay with dead zone.

Figure 26 .
Figure 26.A loci for graphical solution for K = 1.

Figure 28 .
Figure 28.Time responses of closed loop system in Figure 19 for K = 1 with  5 () and relay with dead zone.

Table 2 .
(37)t cycle frequencies for different values of n for Equation(37)where () is an ideal relay.
Figure 14.A Loci for different values of n for Equation

Table 3 .
The frequencies and corresponding amplitudes at arg[()] = −180  as a function of . (

Table 4 .
Recorded values of  and  for K = 1.can then be found by plotting the values given in Table4in (, ) plane as shown in Figure27.From Figure27, it was found that there are two solutions one is (, ) = (3.656, 0.7177 rad s ⁄ ) and the other is (, ) = (0.4308 , 0.2766 rad s ⁄ ).

Table 5 .
The frequency and pulse width of the stable and unstable limit cycles for different values of K.

Table 6 .
Results from simulation, DF and A Locus methods.
A LocusNo solution.System is stable.