Analytical and experimental analyses of nonlinear vibrations in a rotary inverted pendulum

Gaining insight into possible vibratory responses of dynamical systems around their stable equilibria is an essential step, which must be taken before their design and application. The results of such a study can significantly help prevent instability in closed-loop stabilized systems by avoiding the excitation of the system in the neighborhood of its resonance. This paper investigates nonlinear oscillations of a rotary inverted pendulum (RIP) with a full-state feedback controller. Lagrange’s equations are employed to derive an accurate 2-DoF mathematical model, whose parameter values are extracted by both the measurement and 3D modeling of the real system components. Although the governing equations of a 2-DoF nonlinear system are difficult to solve, performing an analytical solution is of great importance, mostly because, compared to the numerical solution, the analytical solution can function as an accurate pattern. Additionally, the analytical solution is generally more appealing to engineers because its computational costs are less than those of the numerical solution. In this study, the perturbative method of multiple scales is used to obtain an analytical solution to the coupled nonlinear motion equations of the closed-loop system. Moreover, the parameters of the controller are determined using the results of this solution. The findings reveal the existence of hardening- and softening-type resonances at the first and second vibrational modes, respectively. This leads to a wide frequency range with moderately large-amplitude vibrations, which must be avoided when adjusting a time-varying set-point for the system. The analytical results of the nonlinear vibration of the RIP are verified by experimental measurements, and a very good agreement is observed between the results of both approaches.


Introduction
The stabilization and tracking of Inverted Pendulums (IPs) undergoing either linear or rotary base motions is a classical engineering problem capable of being employed as a benchmark for various control theories [1]. In this light, the underactuated control [2,3] of these systems, exhibiting inherent nonlinearity and instability, is a challenging subject which, due to its growing application in the development of walking robots, has become the focus of much research during the recent years. Basically, the most common types of IPs are the single-arm Rotary Inverted Pendulum (RIP) [4], the double IP [5], and the cart IP [6,7].
Because of the wide range of applications of IPs, the development of appropriate strategies and methods to control their response has become more appealing and many controllers have been, and still continue to be, tested and introduced in the last few decades. Among these controllers, the linear types are more popular because of their less computational demand to produce control signals, as well as their lower complexity. For example, Proportional-Integral-Derivative (PID) controllers, optimal controllers provided by Linear Quadratic Regulator (LQR) method, and Model Predictive Controllers (MPCs) have been successfully tested and evaluated in some studies [8,9]. In addition, many other well-known linear control algorithms have been broadly tested on IPs by researchers to satisfy required control performances [10,11]. To address the substantial challenges regarding the communication between the controller and the plant, controllers with no bandwidth limitation in comparison with the ones using threshold-based communication have been investigated by Casanova et al. [12]. Also, Hamza et al. have provided a thorough real-time comparison of classical and stateof-the-art control techniques [13]. There also exist more complicated nonlinear controllers capable of global stabilization of IPs. These controllers employ gain scheduling, integral control, feedback linearization, Lyapunov redesign, backstepping, passivelybased control, and high-gain observers to have better performance. Many researchers have attempted to stabilize IPs with the variation and combination of the aforementioned nonlinear controllers, such as Interconnection and Damping Assignment Passively-Based Control (IDA-PBC) [14], cascade non-linear control [15][16][17], Sliding Mode Variable Structure Control (SMVSC) [18], hybrid control for global stabilization [19,20], and nonlinear Model Predictive Controller (MPC) [21,22], to name but a few [23][24][25][26]. However, having an analytical solution-a real-time, low-latency computation-of a closed-loop system offers a clear insight into how variables and control gains affect the result. In addition, since the result of the analytical approximate solution provides insight into the overall behavior of the system, the parameters of numerical procedures-such as time scale-can be determined efficiently. In this regard, quite a few methods have been proposed to find the analytical solution of nonlinear differential equations such as the Hamilton approach [27], the parameter expansion method [28], homotropy perturbation method [29], and other well-known methods [30][31][32]. The exact analytical solution of the nonlinear equations of a large amplitude simple open-loop system pendulum has been found by Sobamowo [33]. The dynamical response of an open-loop 3-DoF pendulum that is harmonically excited has been studied, and the analytical solution implementing the multiple scale method has been carried out [34]. Employing perturbative techniques, Ashrafiuon et al. [35] have presented an approximate analytical solution for the rotary inverted pendulum based on which the step response of the system using the sliding mode controller has been investigated. However, no investigation in the frequency domain has been done yet. In this paper, the analytical solution in the frequency domain has been calculated.
The dynamics of a nonlinear system are significantly richer than that of its linear counterparts, and there are several phenomena that can take place only in the presence of nonlinearity. Research has demonstrated that in the case of an IP the complexity would be even further increased if an elastic beam is used for suspension [36]. Arising from this inherent nonlinearity, the large amplitude vibration of IPs being linearly controlled and locally stabilized in the vicinity of their equilibrium position may lead to instabilities. The linearization of system equations in the neighborhood of the upright position of the pendulum within the design process of the linear controllers restricts their application to the cases with small amplitude disturbances and initial conditions. It necessitates adequate attention to be paid when the set point of the base position is continuously being changed in order to track a predefined path. That is because the existence of a frequency near the natural frequencies of the system in the actuation waveform may result in system instability through the occurrence of resonant large amplitude responses. As a result, although linear experimental model analysis is broadly utilized in the complex structure of many types of inverted pendulums, leading to acceptable results, the substantial increase in computational power and the demand for higher performance of the inverted pendulum is encouraging researchers to take nonlinear effects into account and include nonlinearities in the design and analysis of complex engineering structures.
In the present study, the nonlinear vibratory response of a full-state feedback-stabilized RIP, first introduced by Furuta et al. [37], is both analytically and experimentally investigated. First, using Lagrange's equations, the nonlinear coupled equations of the system incorporating the dynamics of the pendulum and the transfer function of the controller are derived. To verify the later obtained analytical results, a rotary inverted pendulum is fabricated and its geometrical and inertial parameters together with the electrical parameters of the employed DC motor are measured accurately and introduced to the extracted equations. The nonlinear equations are then solved using the method of multiple scales for the case of forced oscillations. The primary resonant response of the system in the neighborhood of both natural frequencies is investigated, validated, and discussed. In addition, controller parameters are set in the process of deriving the solution, and their effect on the frequency response of the system is studied.
2 Analytical investigation Figure 1 shows the schematic view of a rotary inverted pendulum using two rotary encoders as sensors and one DC motor as the actuator. The moving parts of the system may be divided into two groups, where the parts in each group are fixed to each other and may be considered a single rigid body. The first group includes beam 1, rotating parts of encoder 1 (output shaft and the code disk), the rotor of the DC motor, the fixed parts of encoder 2 (the housing and electronics), and all coupling and connections. The second group is composed of the hanging beam (beam 2), the rotating parts of encoder 2, and the couplings used to connect them. Figure 2 clearly shows these two rigid bodies and the parts of the system from which they are composed. This figure includes a 3D model of the fabricated inverted pendulum used for experimental verification of the results which will be introduced in the next section. The total masses of the rigid bodies are represented by m 1 , and m 2 , respectively. In addition, r 1 and r 2 indicate the location of the center of gravity of mass m 2 , where, as shown in Fig. 1, the former denotes the distance of CG from the rotating axis of the motor and the latter is its distance from the rotating axis of encoder 2. To extract the governing dynamical equations of the system, Lagrange's method is adopted. According to this method, first, the total kinetic and potential energies of the system should be determined. The velocity of the center of gravity of m 2 may be written as: The total angular velocity vector of the second rigid body written with respect to the coordinate system xyz is: Therefore, the total kinetic energy of the system will be: where I / designates the moment of inertia of m 1 about the motor axis. Also, I 2x , I 2y , I 2z , and I 2yz are the nonzero elements of the inertia tensor of m 2 about its  Fig. 1 Schematic of the fabricated system mass center and with respect to a coordinate system connected to beam 2, as shown in Fig. 1.
The total potential energy of the system may be written as: By definition, the Lagrangian of the system is: Choosing / and h as our generalized coordinates, the Lagrange's equations of the system will be of the following form: The term Q a in the above equation is the generalized nonconservative force and may be derived by calculation of the total virtual work done: where s denotes the generated mechanical torque of the motor.
For a DC motor, the torque applied to the rotor may be written as a function of its terminal voltage V and speed _ /: where R m and k e are speed regulation and back EMF constants of the motor, respectively. It should be noted that the above equation neglects the inductance of the motor coil which, according to the low frequencies involved in most mechanical systems, is a reasonable assumption.
To make the pendulum stable at h ¼ 0, a closedloop control scheme is used, where the controller accepts /, h, _ /, and _ h as inputs and applies a voltage equal to f c ð/; h; _ /; _ hÞ to the motor. Considering this terminal voltage for the motor and substituting Eqs. (1)-(5), (7), and (8) into (6), the governing equations of the system are derived as: where: The equilibrium points of the system may be found by replacing all terms containing time derivatives of the variable by zero which results in: Therefore, h ¼ 0 and / ¼ / d is an equilibrium point of the system if and only if f c / d ; 0; 0; 0 ð Þ¼0.
As a consequence, assuming the function f c to be linear, it may be written as: It should be noted that, in practice, there may be an error or an offset in the h measurement and, hence, the point at which the measured h vanishes may not exactly coincide with the angle at which CG is accurately located above the axis of the second encoder. Therefore, Eq. (11-a) will not hold anymore, and consequently, h ¼ 0 will not be an equilibrium point for the closed-loop system. To overcome this problem, even if the position angle of the first beam i.e., / is not intended to be controlled at any desired position / d , the term k 1 / must exist in f c . By creating a small error between the angle / and its desired value / d (or with 0 when / d is not important), this term makes Eq. (11-a) hold. Additionally, by employing this term, the system will be able to suspend the pendulum stably at any arbitrary given / d .
Substituting Eq. (12) into (9) while replacing sinh and cosh by their Taylor series around h ¼ 0 and keeping the nonlinear terms of up to third-order results in: The above set of nonlinear equations may be analytically solved using the perturbation method of multiple scales. To do this, assuming the amplitude of angular displacements of / and h to be of order of a small parameter denoted by , one can write: Furthermore, assuming the damping and excitation terms to be of order of 3 , one can let: New independent variables called time scales may be introduced accordingly as: Representing the operator denoting partial differentiation with respect to T n by D n (i.e., D n ¼ o oT n ) and using the chain rule of differentiation, one can write operators denoting first and second derivatives with respect to t as an expansion of D n s and their multiplications as: Substituting Eqs. (14), (15), and (17) into (13), then separating terms of different orders of , one has: e 3 : where The general solution of Eq. (18) may be written as: where cc denotes the complex conjugate terms. The above solution assumes that the following equation from which x 1 and x 2 may be calculated, has four real roots: This condition needs k 1 and k 2 to be chosen so that the following inequalities hold: The coefficients a 1 and a 2 in Eq. (21) are as follows: To analyze resonance occurrence in the neighborhood of the natural frequencies of the system, one may let: where X 1 and X 2 differ from x 1 and x 2 by small parameters r 1 and r 2 , respectively, as: To eliminate the secular terms in Eq. (18), the coefficients of the terms expðix i T 0 Þ in the following expressions must vanish: where Substituting Eqs. (21), (25), and (26) into (27), one obtains: where Eliminating the excitation and nonlinear terms from the above expressions, the condition for linear stability of the system is derived as: Amplitudes A 1 and A 2 may be represented by polar form as: Substituting the above equation into Eq. (29), one has: By the separation of real and imaginary parts in Eq. (33) and introduction of the following new variables: the following set of autonomous equations are found: The equilibrium points of the above dynamical system may be found by substituting a The above equations may be solved to calculate vibration amplitudes a 1 and a 2 as functions of the excitations amplitudes E 1 , E 2 ; and deviations r 1 and r 2 from the natural frequencies of the system. To depict the response curves, it suffices to simply calculate r 1 , r 2 , E 1 , and E 2 by a rearrangement of Eq. (36) as: To explore how the inherent geometric nonlinearity of the pendulum affects the behavior of the system, a linear solution to the governing equations is obtained as well. The results of the linear solution are used along with the results of the nonlinear solution when comparing the analytical, numerical, and experimental findings in the next section. The first step toward this linear solution is eliminating the nonlinear terms in Eq. (13) which results in: where / d is a single harmonic excitation as: Substituting the above expression into Eq. (37), the steady-state solutions for / and h are found to be of the following form: where U and H are complex amplitudes, derived from the following equation: Solving the above equation, one obtains:

Experimental setup
A photo of the fabricated rotary inverted pendulum with the aim of verification of the analytical results is shown in Fig. 3. The measured and calculated values for the parameters of the system are as given in Table 1. The electric parameters of the motor are determined by measuring the input current in both noload and locked-rotor conditions. To obtain an accurate estimation of the structural parameters of the system, they are calculated using Solidworks. To this end, first, a precise 3D model of the system is generated with corresponding materials assigned to all the parts, and then the total mass, moment of inertia, and the location of center of mass is calculated for the parts forming rigid bodies 1 and 2 (See Fig. 2). All the analytical and numerical results presented in the next section are for a system with parameters given in Table 1.

Results verification and discussion
In this section, we present and compare the analytical, numerical, and experimental results. At first, the acceptable values of k 1 and k 2 for which the system is stable are determined based on inequalities of Eq. (23). By depicting the corresponding curves, one may find the acceptable area illustrated in Fig. 4. k 1 and k 2 can also be found either experimentally or through any other methods of control such as poleplacement or LQR. After finding an appropriate set of values for k 1 and k 2 to be 50 and 250, respectively, the natural frequencies of the present 2-DoF system can be calculated by substituting the obtained values of k 1 and k 2 along with the measured parameters listed in Table 1 into Eq. (22). By doing so, the natural frequencies of the system are found to be x 1 ¼ 4:24 rad=sec and x 2 ¼ 50:2 rad=sec. Figure 5 illustrates the acceptable area for k 3 and k 4 found by depicting the inequalities of Eq. .
Acceptable values of these state feedback gain led to an asymptotically stable system. Figures 6, 7, 8, 9 show the analytically-calculated frequency responses of the system including both angular displacements h and / in the neighborhood of the first and second resonances for different amplitudes of excitation. To verify the performed analytical solution, the governing equations are also numerically solved employing the Runge-Kutta method. The results of this numerical solution are depicted together with the analytical resonance curves in Figs. 6, 7, 8, 9 for comparison purposes, where a very good agreement is observed. As these figures demonstrate, at small amplitudes of excitation the resonance curves are almost symmetric, resembling the behavior of a linear oscillator. As excitation amplitude increases, this symmetry gradually disappears through bending of the curves. As shown by the results, the nonlinearity has a hardening effect on the response curve in the neighborhood of the first resonance, while at the second resonance the system shows a softening behavior and the curve is bended toward lower frequencies. Another important point which can be inferred from Figs. 6,7,8,9 are the different levels of agreement between analytical and numerical results in different frequency bands. Although these two solutions are in good agreement in the vicinity of the resonance peak, they start to diverge as the excitation frequency moves away from the natural one. This occurs because the perturbation solution assumes the detuning parameter r to be of order 2 (See Eq. 26) and  therefore is valid only in a close neighborhood of the natural frequency. One of the significant characteristics of these nonlinear frequency responses is the existence of a multivalued region which emerges and widens as the excitation amplitude increases. In this region, the curve is composed of three branches including one middle branch with saddle stability enveloped by two low and high energy stable branches. There are two bifurcation points of the saddle-node type in which a saddle point and a stable node of the system collide and annihilate each other. The location of these bifurcation points indicates the edges of the multivalued interval and the trajectory of the excitation (trajectory of / d in this case) determines whether the response falls on the upper branch of the curve or the lower one in this region. Figure 10 clearly demonstrates the evolution of a fully nonlinear response curve from an almost symmetric linear-like one as the excitation amplitude continuously grows.
To demonstrate the interaction of vibration modes corresponding to different natural frequencies, vibration amplitudes are depicted as surfaces that are functions of both r 1 and r 2 in Figs. 11 and 12. As  these figures illustrate, while a 2 is affected by frequencies of both excitation terms (See Eq. 25), the variation of a 1 as the frequency of second mode excitation changes is negligible. The effect of excitation amplitudes on the amplitude of the response at both modes of vibration is shown in Figs. 13 and 14. Additionally, Fig. 15 through Fig. 18 show the amplitude of the response for both h and / with respect to the amplitude of excitations for several values of r 1 and r 2 which, as expected, contain both multivaluedness and jump phenomenon, resulting from the nonlinear phase-amplitude interaction governed by Eq. (35).
To verify our analytical and numerical solutions, the frequency response curves of the system for both h and / at the first resonance are found through experimental measurements as well. The results are provided in Figs. 19 and 20 together with the numerical and analytical solutions. It is worth noting that the experiments are conducted at constant amplitude of excitation, while the angular frequency is increasing with the step size of 0:2rad=s. The analytical linear Vibration amplitude a 1 with respect to r 1 and E 1 (state feedback gains: k 1 ¼ 50, k 2 ¼ 250, k 3 ¼ 4, k 4 ¼ 12) Fig. 11 Vibration amplitude of a 2 with respect to r 1 and r 2 (state feedback gains: k 1 ¼ 50, k 2 ¼ 250, k 3 ¼ 0:8, k 4 ¼ 0:8; Excitation amplitudes: E 1 ¼ 0:03 and E 2 ¼ 0:03) solution as derived in Eq. (42) is depicted together with the nonlinear response in these figures to demonstrate how the experimental results deviate from the symmetric curve resulting from the performed linear analysis. As can be clearly seen, in the close neighborhood of the natural frequency the linear solution leads to a considerable error, while the nonlinear one resembles the experimental behavior to a large extent.

Conclusion
In this paper, the nonlinear vibration behavior of a linearly controlled rotary inverted pendulum (RIP) is investigated through analytical simulations and Fig. 12 Vibration amplitude of a 1 with respect to r 1 and r 2 (state feedback gains: k 1 ¼ 50, k 2 ¼ 250, k 3 ¼ 0:8, k 4 ¼ 0:8; Excitation amplitudes: E 1 ¼ 0:03 and E 2 ¼ 0:03) Fig. 13 Vibration amplitude of a 1 with respect to E 1 and E 2 (state feedback gains: k 1 ¼ 50, k 2 ¼ 250, k 3 ¼ 0:5, k 4 ¼ 0:7; Deviations from Natural frequencies: r 1 ¼ 0, r 2 ¼ 0) Fig. 14 Vibration amplitude of a 2 with respect to E 1 and E 2 (state feedback gains: k 1 ¼ 50, k 2 ¼ 250, k 3 ¼ 0:5, k 4 ¼ 0:7; Deviations from Natural frequencies: r 1 ¼ 0, r 2 ¼ 0) Fig. 15 Amplitude of h with respect to E 2 with several values of r 2 (state feedback gains: k 1 ¼ 50, k 2 ¼ 250, k 3 ¼ 0:03, k 4 ¼ 0:25) Fig. 16 Amplitude of h with respect E 1 with several values of r 1 (state feedback gains: k 1 ¼ 50, k 2 ¼ 250, k 3 ¼ 0:03, k 4 ¼ 0: 25) experimental measurements. The perturbation method of multiple scales is used to solve the governing coupled nonlinear equations of the system derived using Lagrange's variational principle. To experimentally verify the predicted results of the analytical model, both of its mechanical and electrical parameters are set to be equal to those of the real RIP system employed in the experiment. The comparison of the results obtained through the analytical, numerical, and experimental approaches reveals a good agreement among them. The acceptable range of values for the controller feedback gains, leading to the stabilized oscillatory motion of the pendulum, is derived through the performed analyses. The findings indicate that, while the nonlinearity involved in the first mode of the vibration is a hardening type, the system shows a softening behavior when excited in the vicinity of its second natural frequency. As the amplitude of the alternating part of the desired position waveform is enhanced due to the bending of the resonance curves, large oscillations are possible to occur within a wider frequency band around the natural frequencies. This calculated frequency band must be avoided in practice when the set-point of the system is being changed as time passes, particularly in the case of periodic motions.