Two-Parameter Dynamics of An Autonomous Mechanical Governor System With Time Delay


 Understanding of dynamical behavior in the parameter-state space plays a vital role in the optimal design and motion control of mechanical governor systems. By combining the GPU parallel computing technique with two determinate indicators, namely, the Lyapunov exponents and Poincaré section, this paper presents a detailed study on the two-parameter dynamics of a mechanical governor system with different time delays. By identifying different system responses in two-parameter plane, it is shown that the complexity of evolutionary process can increase significantly with the increase of time delay. The path-following strategy and the time domain collocation method are used to explore the details of the evolutionary process. An interesting phenomenon is found in the dynamical behavior of the delayed governor system, which can cause the inconsistency between the number of intersection points of certain periodic response on Poincaré section and the actual period characteristic. For example, the commonly exhibited period-1 orbit may have two or more intersection points on the Poincaré section instead of one point. Furthermore, the variations of basins of attraction are also discussed in the plane of initial history conditions to demonstrate the observed multistability phenomena and chaotic transitions.


Introduction
Time-delay effect is widely encountered in many dynamical systems including biological systems [1][2][3], social and economic systems [4,5], neural networks [6,7], chemical reactions [8,9], and mechanical systems [10][11], due to the essential dependency of present states on past states. Compared to the dynamical systems without time-delay effect, systems involving inherent time delays can exhibit more rich and complicated dynamics, including the instability of equilibrium, multiple stability switches induced by time-delay [12,13], high co-dimension bifurcations [14][15][16], coexistence of attractors [17,18], and chaotic motions [19,20]. In addition to the understanding of delay-induced bifurcations and transitions, the beneficial applications of time-delay effect have also attracted great interest from the research community. For example, time-delay effect has been considered in the design of real-world controllers, due to the unavoidable time consumption in the measurement, signal processing, computing and transmission processes. The inherent time delays could be used to achieve the desired control objectives, such as for the traffic arrangement [21,22], signal enhancement [23], broadband energy harvesting [24,25], and absorber and isolator design [26,27].
The mechanical centrifugal governor has played a role of milestone in the field of automatic feedback control, since it was invented and used for stabilizing the working speed of steam engines. Owing to the simple yet effective characteristics, the centrifugal governor devices have been widely used in many rotating machines to avoid potential damage, such as in the internal combustion engines, fueled turbines, modern striking clocks and so on. Induced by the inherent nonlinearities, many dynamical phenomena have been observed and reported in such systems, including periodic oscillations [28,29], chaotic motions [30,31], hyperchaotic behavior [32], multistability [15,33], fractal and mode-locking structures [34,35]. To avoid the damages caused by the undesired system response, various control schemes have been proposed to stabilize the system. For example, the feedback controller [36,37], adaptive controller [37], and time-delayed controller [38] were used to convert the chaotic motions of the governor systems into specified motions. A synchronization-based adaptive sliding mode controller [39] and a robust adaptive controller [40] were designed to realize the finite-time stabilization of mechanical governor systems with uncertainties. The design of anti-bifurcation controllers [41,42] were studied in a centrifugal governor system to generate desirable oscillations which can be applied in the fault diagnosis, monitoring, mixing and low-energy navigation.
Although remarkable achievements have been made to the understanding of dynamical behavior in the literature, few research works have been devoted to the dynamics of such governor systems under the inherent effect of time-delay. Actually, time-delay can be originated from the gear transmission, spring deformation [37], steam transportation or chemical combustion [15]. The reaction delay caused by the mass-spring structure was considered in the governor system, and varieties of dynamical responses were illustrated from numerical simulations [37]. In Ref. [15], the time-delay effect in the governor system was modeled by considering the reaction delay induced by the working principle of engines. It was found that the time delay can shrink the stability region of equilibrium and even induce high co-dimension bifurcations, e.g., resonant and non-resonant double Hopf bifurcations. However, the dynamics and transitions are not fully understood and yet to be studied in the parameter region, where the chosen equilibrium loses its stability. To develop an in-depth insight into the open questions unaddressed in Ref. [15], this paper aims to reveal the global dynamics of the autonomous governor system with and without time delay, e.g., the two-parameter classification of system response and the related multistability phenomena, which will develop a deep understanding of time-delay induced global dynamical behavior of the present governor system.
Indeed, the general studies on the nonlinear dynamics were mainly focused on the local bifurcation analysis during past decades, which only provides a limited view on the system performance. With the development of computer technology, it has been possible for the researchers to study a wide range of bifurcation behaviors of certain systems in the parameter and state spaces. It is worth mentioning that an efficient accelerated tool, i.e., the CUDA-based GPU parallel computing technique [43], has been rapidly developed for the computation of multi-parameter dynamics of the systems and the efficiency of such a numerical tool has been proved in many works [44][45][46]. However, these works are basically directed to the dynamical systems described by the ordinary differential equations, few research works are associated with the global dynamics computation for the dynamical systems governed by the delayed differential equations (i.e., functional differential equations). In this paper, the CUDA-based GPU parallel computing technique is adopted to study the classification of two-parameter dynamics in a time-delayed governor system. Therefore, this work can be regarded as an important attempt in the computation of global dynamics of the delayed dynamical system.
Given that the existing numerical methods are basically for the ordinary differential dynamical systems, the treatment of the delayed terms and transforming the delayed differential dynamical system into the ordinary differential dynamical system then become a noteworthy concern in the computation of global dynamics of the delayed systems. Currently, the efficient approximation schemes for the delayed dynamical systems include the Galerkin projections [47] and the discretization approaches [48]. Based on these approximation schemes, the delayed dynamical system can be transformed into a high-dimensional ordinary differential system. Significant time-consumption will be involved in precisely obtaining the portfolios of global dynamics of such a high-dimensional system in the parameter and state spaces.
In addition to the high time-cost on the computation of global dynamics, the delay-induced winding phenomenon in the identification of bifurcation will be another challenge. As reported in Ref. [49], the self-intersection phenomenon of phase trajectory can be delay-induced in the autonomous and delayed dynamical system, which will lead to the difficulties in distinguishing periodic response and identifying the bifurcation mechanism.
For the present autonomous and delayed governor system, the winding phenomena of periodic response are observed with certain delays and the bifurcation identification becomes a difficult task. To address this problem, a type of semi-analytical method, i.e., the time domain collocation method (TDCM) [50], is adopted to determine the bifurcation types of periodic response for the following reasons. First, the TDCM is easily implemented without the redundant process of formula derivation and the limitation of nonlinear terms of dynamical systems.
Second, to the best of the authors' knowledge, this method has only been applied in the analysis of the low-dimensional dynamical systems with simple nonlinearity [50][51][52]. It is expected that this method will also be effective in solving the high-dimensional dynamical system, especially for the delayed dynamical system with complex nonlinearity. Based on these two reasons, the TDCM is chosen in this study to obtain the approximate analytical solution of periodic response and identify corresponding bifurcations.
The main contributions of this work can be briefly summarized in three aspects.
Firstly, for the mechanical governor system, it is necessary to study the global dynamics of the system in the parameter and state spaces with inherent time delay for a better understanding of the system performance. Secondly, for the overall development of the delayed nonlinear dynamical systems, the exploration on the computation and bifurcation identification of global dynamical behaviors in this paper will provide a meaningful reference to the numerical simulation of complicated dynamical behaviour. Thirdly, the winding phenomena for periodic responses are observed and further discussed in the delayed governor system. It is found that the Poincaré section analysis may fail its effectiveness in distinguishing period property of periodic orbit in some cases.

Governing equations of the mechanical governor system
As developed in Ref. [15], the governing equations of a type of mechanical governor system with reaction time-delay generated from the working principle of engines can be described by the following functional differential equations: 22 ( ) sin cos cos sin cos dx y dt dy e pz and the detailed modelling process can be found in Ref. [15]. For the delayed governor system (1), one equilibrium can be easily determined as: By adopting the coordinate transformations can be re-expressed as For the mechanical governor system (3), the linear stability and local bifurcation behavior as well as the effect of time delay have been well investigated in a two-parameter plane, i.e., the ( , ) bQ plane in Ref. [15]. As an extension of Ref. [15], this paper will investigate the global behavior of the time-delayed system and discuss the influence of time delay in the same two-parameter plane. More specifically, the system parameters are identical to Ref. [ can be found in Ref. [15].
For the sake of description, rewriting the governing Eq. (3) into a compact form

Numerical integration algorithm
For the dynamical system with time delay and strong nonlinearity, the exact analytical solution is almost unavailable and the direct numerical integration becomes an indispensable technique to solve the complex functional differential equations. In this study, the fourth-order Runge-Kutta method with fixed step and the linear interpolation technique are combined to solve the delayed differential Eq. (4) with the initial value problem.
To begin with, the procedure of dimension reduction of the delayed governor system (4) is implemented via the discretization operation [48]. By dividing the initial , where N  is a large positive integer. Based on the discretization operation, the continuous initial function () t  can be approximated as an initial vector 0 Y : Let h be the integral step size and satisfy the condition ht  , the single-step iteration map  can then be conducted as: By selecting the proper numerical integration algorithm, all components of the new iteration vector h Y can be expressed as: is obtained by using the linear interpolation,

Largest Lyapunov exponent and Poincaré section
Lyapunov exponent is an effective indicator for the quantitative classification of different system responses. For an autonomous dynamical system, the system is in the state of equilibrium if the largest Lyapunov exponent (LLE) is negative. Chaotic oscillation will appear in the system response if there exists a positive LLE. Otherwise, the system will exhibit periodic or quasi-periodic motions, if the LLE equals to zero [53]. When the periodic and quasi-periodic motions need to be further distinguished, the second largest Lyapunov exponent or the Poincaré section will be examined. The LLE has received the most attention for its high-efficiency in distinguishing different dynamical responses. In this paper, the estimation of LLE is based on the method proposed in Ref. [54]. According to the estimation scheme [54] as shown in Fig. 1, the largest Lyapunov exponent is calculated by where h is the integral step size and 1 n is the total number of iterations. Symbol 0 d represents the initial distance between two nearby phase trajectories, i.e., (0) 0 Y and (0) 1 Y . i d represents the evolution distance between the two trajectories at the i th iteration. In this scheme, renormalization is required for each iteration, namely In addition to the Lyapunov exponents, the Poincaré section is another widely used tool to distinguish the types of system responses. Generally, the system exhibits chaotic oscillation when a large number of disordered discrete points are found on the Poincaré section. The system will show quasi-periodic oscillation when a closed curve is observed on the Poincaré section. Periodic oscillation occurs in the system response when limited and countable points are observed on the section. However, it might be inaccurate to determine the period property of periodic response by only counting the points on the Poincaré section for the autonomous system with time delay, since the self-intersection phenomenon of phase trajectory [49] may mislead the classification.
To overcome this confusion, the other methods will also be applied and described in subsequent anlaysis. In this paper, the Poincaré section is selected as the section of zero angular velocity and given by: It should be mentioned that the point of zero angular velocity on Poincaré section is detected from the binary process by varying integral step size h . Moreover, the Poincaré section  is mainly used for the classification of different periodic responses in two-parameter diagram and the presentation of single-parameter bifurcation diagram.

GPU parallel computing
Using the above mentioned numerical computation schemes, the single-parameter bifurcation diagram of the delayed mechanical governor system can be easily obtained. However, it is often a partial understanding of the delayed governor dynamics by only varying a single parameter. To explore the variation of the dynamical behavior of delayed dynamics with simultaneous changes of multiple system parameters, the GPU parallel computing technique is considered here to achieve the two-parameter diagram.
By dividing the two-parameter plane ( , ) bQ into LL  grid points uniformly,

Time domain collocation method
In the study, the time domain collocation method (TDCM) [50] is chosen as a theoretical tool to identify the bifurcation types of periodic response. The principle of TDCM will be briefly described in this subsection. Since the mechanical governor system (4) is autonomous and the frequency is unknown, the time re-scaling operation of Eq. (4) is conducted firstly for the convenience of the use of TDCM. Let 1 , tt      , then the time re-scaling governor system becomes:  . The period of system (11) has been transformed into 2 at this time. Approximating the solution 1 the Fourier series expansion yields: coefficients and N is the harmonic number in the truncated Fourier series solution.
Correspondingly, the derivative and the delayed state can be expressed as: Substituting Eq. (12)-Eq. (14) into the re-scaled governing Eq. (11) leads to the expression of the residual function: Thereafter, a series of residual functions for all the equidistant points are derived as: To obtain the above 3 (2 1)+1 N  unknowns, i.e., 3

Transition matrix for stability analysis
To analyze the stability and bifurcation of the detected periodic response, the state transition matrix will be needed. Based on the obtained harmonic solution of system (11) and regarding it as the initial condition, the state transition matrix i G in single integral step can be theoretically determined as where N  represents the segments derived from the equidistant division of the initial   Figs. 2(a-f) with those in Ref. [15] shows that the stable regions of equilibrium and its boundaries are identical to each other.
Additionally, it has been reported in Ref. [15] that the equilibrium would lose its stability mainly through Hopf bifurcation. The bifurcating limit cycle (indicated by symbol N1 and pink regions) are also identified near the stability boundaries of the equilibrium states, as shown in Figs. 2(a-f). Moreover, the fine structures lying on the left bottom of the two-parameter plane and the steady-state response originating from a certain co-dimension two bifurcation of equilibrium, which have been reported in Ref. [15], are also detected and shown in Figs. 2(a-f). These features indicate the effectiveness of the numerical computation strategies designed in this paper.
For the nonlinear dynamics uncovered in Ref. [15], the numerical results

presented in Figs. 2(a-f) can now provide a clear interpretation and bring new insights.
For example, the N2 N6 : types of periodic responses and chaotic motions revealed here have not been reported in Ref. [15]. In the parameter region of chaotic motions in  Fig. 2(a), the single-parameter bifurcation diagrams are obtained by using the direct numerical integration (DNI) and path-following strategy.
The corresponding results in terms of the chosen Poincaré section determined by Eq.  Fig. 4(a), the appearance and disappearance of N3 type of periodic response are both induced by the SN bifurcation. However, some coexisting attractors can be newly discovered and displayed in the inserted windows in Fig. 4(b). These newly detected attractors might be often neglected if the continuation analysis of N3 type of periodic response is not fully considered. By increasing the bifurcation parameter b of Fig. 4(b), the stable N3 type of periodic response shall lose its stability through the subcritical period-doubling (PD) bifurcation on the red dotted line, as shown in the inserted windows. The unstable N3 type of periodic solution marked by the red dots appears accordingly. Thereafter, the unstable N3 type of periodic solution will regain its stability by a supercritical PD bifurcation. What should be noted is that chaotic attractor will appear near such a supercritical PD bifurcation. And the newly stable N3 type of periodic response will finally disappear via an SN bifurcation.
For the case of 1   of Fig. 2(b), the single-parameter bifurcation diagrams and the continuation results are presented in Figs. 5(a-b). From Fig. 5(b), it can be seen that a typical PD bifurcation route from periodic solution to chaotic solution is observed. Therefore, the focus will be put on the jump phenomena in Fig. 5(a).
Different from the jump mechanism of N1 type of periodic response revealed in Fig.   4(a), the case in Fig. 5(a) is more complicated. Two branches of unstable periodic solutions, i.e., branches 1 and 2, are traced. To better illustrate the evolutionary process in the case, variation tendencies of the largest two Floquet multipliers with the change of parameter Q are presented in Fig. 5(c-d). Here, the trivial Floquet multiplier 1 (which always exists) is not displayed in the diagrams. Apart from this trivial multiplier 1, the rest multipliers are ranked and the top two multipliers, i.e., 12 , , are displayed. The ranked multipliers corresponding to the branch 1 are shown in Fig. 5(c), while those of branch 2 are shown in Fig. 5(d). It can be seen that a limit point exists at the intersection point between these two branches, since a nontrivial Floquet multiplier 1  is calculated and the largest Floquet multiplier is larger than 1, i.e., 12 Fig. 5(d), which indicates a NS bifurcation of branch 2.
As for the jump between N3 type of periodic motion and chaos in Fig. 5(a), similar scenarios to Fig. 4(b) are found and the newly bifurcating system responses are shown in the inserted windows.
From the above analyses, rich phenomena of multistability are observed in the system response, such as the coexistence of periodic response and chaotic motion, and coexistence of varieties of periodic responses. The details on these coexisting response will be discussed later.

Winding phenomena and bifurcations with 2  
When increasing the time delay to 2   , the single-parameter bifurcation diagrams are plotted in Figs. 6(a-b). It is easy to notice that a supercritical PD bifurcation will arise in these cases. And the system will exhibit simple behaviors, e.g., the equilibrium state and simple periodic oscillations. When the time delay is further To address this problem, the TDCM will be applied to obtain the analytical solutions of the complex periodic orbits and further clarify the bifurcation types.
Taking region I (black dotted circle) in Fig. 7(a) as an example, 12 N  is chosen in the Fourier series expansion, i.e., Eq. (12), and the number of collocation points M is set as 34 M  in this case. Accordingly, the analytical solutions of the periodic responses corresponding to region I can be derived and marked as red circles in Fig. 7(c). The relevant largest two Floquet multipliers are exhibited in Fig. 7(d).
From Fig. 7(c), it can be seen that the theoretical results are in a good agreement with the numerical results. It is also noticed that the PD bifurcation occurs on the red dotted line and no bifurcations are detected on the black dotted lines by analyzing the Floquet multipliers in Fig. 7(d). To explore the transition mechanism of just beforeand after-the black dotted lines in Fig. 7(c), the relationship between the bifurcation parameter and the harmonic components of analytical solutions, i.e., the harmonic orders and amplitudes, is studied next.
To begin with, re-expressing the harmonic amplitudes in Eq. (12) as Besides, the fact that the existence of even harmonic components in the analytical solutions can induce the asymmetry of periodic orbits [56] is also witnessed in this analysis. Same as the bifurcation identification procedure for the complex periodic response in Fig. 7(a), the analytical solutions and the Floquet multipliers are also calculated for the case of Fig. 7(b), and a PD bifurcation is determined in this case.  Fig. 8(c). Meanwhile, the coexistence of attractors can be observed in this region, as shown in the inserted phase diagram. For the bifurcations in region Ⅲ , the N1 and N3 types of periodic responses will both disappear by the SN bifurcation, since the approximated Floquent multiplier 1  can be predicted in Fig. 8(d-e). In region Ш, the coexistence of these independent periodic solutions can be observed and the coexisting attractors are shown in the inserted phase portrait of Fig. 8(a). For the bifurcation diagram shown in Fig. 8(b), the supercritical PD bifurcation route arises. However, the evolutionary process of periodic responses shown in Fig. 8(b) is obviously more complicated than the case in Fig. 7(b).
For the case of 5   , the resultant single-parameter bifurcation diagrams are displayed in Figs. 9(a-b). As mentioned previously, the evolutionary process in this case is much more complicated than that in Figs. 8(a-b). That is, the relationship between the distribution of intersection points on the chosen Poincaré sections and periodic characteristic becomes more indistinguishable. To distinguish different responses and clarify bifurcation behaviors in these cases, the Floquet multipliers are also calculated and the corresponding bifurcation types are inserted in Figs. 9(a-b).
The PD bifurcation is the main bifurcation type in these two cases. Specifically, a supercritical NS bifurcation tree is found in Fig. 9(a), which coexists with a PD bifurcation tree. Accordingly, the multistabilities, i.e., periodic response coexisting with quasi-periodic response, can be observed. Moreover, the coexistence pattern of two independent periodic responses is detected in region Ⅳ of Fig. 9(a) and will be discussed in the following subsection.

Co-dimension two bifurcation of equilibrium for 5  
A co-dimension two Hopf bifurcation point of equilibrium could be found in the neighborhood (in terms of parameters) of region Ⅳ of Fig. 9(a) and has been reported in Ref. [15]. Based on the results of Ref. [15], four types of responses will arise in region Ⅳ, i.e., equilibrium 0 E , two periodic solutions 1 E and 2 E , and one unstable quasi-periodic solution 3 E . The bifurcations and transitions between the four responses are schematically presented in Fig. 9(c). Using the DNI and path-following methods, the steady response of system in region Ⅳ is traced and indicated by blue dots in Fig. 9(d). It can be seen that the transitions of Fig. 9(d) are in accordance with the transitions marked by the black solid lines in Fig. 9(c). As for the jump phenomenon arising between 1 E and 2 E in Fig. 9(c), it is also verified by numerical results in Fig. 9(d).
To verify the transition mechanisms in terms of the periodic solutions 1 E and 2 E , the analytical solutions are determined as the red circles shown in Fig. 9

Multistable phenomena and chaos routes
This subsection discusses the qualitative or quantitative changes of the system responses in the state plane. As mentioned previously, a number of coexisting attractors has been observed in the mechanical governor system with and without time delay. For the non-delay cases shown in Figs. 4(a-b), two main bifurcation trees are observed under different initial perturbations. One bifurcation tree originates from the supercritical PD bifurcation of N1 type of periodic response and the other is caused by the supercritical PD bifurcation of N3 type of periodic response. Accordingly, rich coexistence phenomena can be induced in the governor system without time delay, such as the coexistence of different periodic attractors, coexistence of different periodic attractors and chaos. In addition to the above-mentioned coexistence patterns, the coexistence of equilibrium and periodic response can also be observed in this case.
To study the evolutionary process of basins of attraction, the coexistence of N3 type of periodic response and chaos arising in Fig. 4(a) is taken as an illustrative example. Four groups of parameters are selected to compute the basins of attraction and the results are presented in Figs. 10(a-d) accordingly. The basin of attraction is also obtained by the GPU parallel computing technique. The state space is considered as follows, 3 1 {( , , ) For the mechanical governor system with time delay, some coexistence phenomena are newly found. When the time delay is small as the case shown in Fig.   5(a), a new coexistence of stable N1 type and N3 type of periodic responses emerges.
As the time delay increases, new dynamic phenomena are observed, e.g., the coexistence of two branches of NS bifurcation tree in region Ⅱ of Fig. 8(a) and the co-dimension two bifurcation of equilibrium in Fig. 9(a). Accordingly, two coexistence patterns, i.e., coexistence of N1 type and N3 type of periodic attractors, coexistence of periodic and quasi-periodic attractors, occur in the cases of Fig. 8(a) and Fig. 9(a), respectively. To illustrate the coexistence of the N1 type and N3 type of periodic attractors for the case of Fig. 8(a), the basin of attraction is considered to be computed in the following sate space: In the case, the state plane 00 ( , ) uw is displayed and also divided into 300 300  grid points. Moreover, four groups of parameters are chosen to study the evolution of the basin of attraction. The corresponding results are exhibited in Fig.   11, where the red domain is the basin of attraction of the N3 type of periodic response, and the green domain represents the basin of attraction of N1 type of periodic response. The rest rufous domain corresponds to the diverging motion under the chosen initial condition. Same as the case in Fig. 10, the self-similar basin structure is also observed in Fig. 11. Due to the range limitation of the selected state plane, the self-similar basin structures located on the right hand sides of Figs. 11(a-d) are partially presented and would not be further studied in this paper. It is easy to find that the evolution tendency of the basin of attraction of the two different periodic solutions is identical with the bifurcation diagram in Fig. 8(a).
In the chaotic transitions exhibited in the above results, three main routes leading to chaos can be found in the system response. The first route is that the chaotic motion is generated from a PD bifurcation and then disappears by the way of boundary crisis [57], e.g., Figs. 4(a-b) and Fig. 5(a). In these cases, the PD bifurcation is involved with the supercritical type. The second route is that the chaotic motion is born and destroyed both via a PD bifurcation, e.g., Fig. 4(b), Fig. 5(b) and Fig. 7(a). The third route is that the chaotic motion arises from a NS bifurcation, such as Fig. 8(a). In Fig.   8(a), the subcritical NS bifurcation of N1 type of periodic response is involved and such a route has ever been reported in Ref. [57].
To explore the global behavior of the governor dynamics with and without time delay, a two-parameter dynamic analysis has been carried out by using the CUDA-capable GPU computing technique. The uncovered dynamics in Ref. [15] has been systematically investigated by the single-parameter bifurcation diagram and path-following technique. It was found that the complex periodic orbits could arise in Based on the work presented in Ref. [15] and the present study, the two-parameter dynamics of the delayed governor system has been extensively investigated. It is anticipated that these research results would provide a meaningful reference to the optimal design and motion control of the mechanical systems and the exploration of local and global dynamics of the delayed systems.

Declarations Funding
This work has been supported by the Key Program of National Natural Science        The top two nontrivial Floquet multipliers for the N3 type periodic response in region Ⅱ of Fig. 8(a), (d-e) The top two nontrivial Floquet multipliers for the N1 and N3 type periodic responses in region Ⅲ of Fig. 8(a), respectively.