Dynamic response and chaotic behavior of a controllable flexible robot

Flexible robots with controllable mechanisms have advantages over common tandem robots in vibration magnitude, residual vibration time, working speed, and efficiency. However, abnormal vibration can sometimes occur, affecting their operation. Traditionally only simple mechanisms are considered in studying abnormal vibration, omitting reference to important chaotic phenomena caused by the flexibility of the mechanism rod. In order to better understand the causes of abnormal vibration, our work takes a controllable flexible robot with a complex series-parallel mechanism as a research object and uses a combination of Lagrangian and finite element methods to establish its nonlinear elastic dynamics. The effectiveness of the model is verified by comparing the calculated frequency with the frequency measured in a test. The time-domain diagram, phase diagram, Poincaré map, maximum Lyapunov exponent, and bifurcation diagram of the elastic motion of the robot wrist are studied, and the chaotic phenomena in the system are identified through the phase diagram, Poincaré map, the maximum Lyapunov exponent, and the bifurcation diagram. The relationship between the parameters of the robot motion and the maximum Lyapunov exponent is discussed, including trajectory angular speed and radius. The results show that chaotic behavior exists in the controllable flexible robot and that trajectory angular speed and radius all have an influence on the chaotic motion. Our work provides a theoretical basis for further research on the control and optimal design of flexible robot mechanisms.


Introduction
Controllable mechanisms include adjustable mechanisms, hybrid drive mechanisms, variable input parameter mechanisms, and flexible robots. Controllable flexible robots have the following advantages over traditional robots: higher intelligence, greater flexibility of output motion, and higher accuracy of motion. They belong to a class of highly coupled nonlinear systems with complex nonlinear dynamic behaviors. Chaos is one of the most important dynamic characteristics of a nonlinear dynamic system. It is a bounded irregular motion distinct from periodic motion and quasiperiodic motion. It can well describe the behavior of nonlinear systems. The elastic deformation of mecha-nism rods can cause the mechanism to behave chaotically. At present, the methods of analyzing chaotic motion of a system include the time history method, the phase diagram method, the Poincaré mapping method, the maximum Lyapunov exponent method, and the bifurcation diagram method [1][2][3][4], etc. It is therefore important to study the dynamic response and the chaotic behavior of controllable flexible robots, and establish the dynamic performance design theory of this type of robot, which will help establishing further characteristics of the nonlinear dynamic behavior of this system and in revealing its abnormalities [5][6][7][8].
In the past few decades, many researchers have worked on the dynamic response of robots. Tang et al. [9], based on Hamilton's principle, established a set of nonlinear equations of motion of a double-cable singlelink flexible manipulator with two boundary constraints and studied the natural frequency and terminal response of the system. The effectiveness of the cable in suppressing the vibration of the system was verified by both simulation and experiment. Cai et al. [10] proposed a multi-objective optimization design method to improve the elastic strain response distribution of a joint strain gauge in the working area, and conducted related experimental studies to improve the sensitivity and resolution of the joint torque sensor, finally achieving compliance control of the robot joints. Amer et al. [11][12][13][14][15] used a Lagrangian method to establish nonlinear dynamic models of single-pendulum, doublependulum, and triple-pendulum systems. A multi-scale method was used to solve the problem analytically. The solvability condition was checked and the time history and resonance response of the solution were obtained. Using a nonlinear stability analysis method, the stable and unstable regions were analyzed, and the influences of different physical parameters on dynamic motion were investigated. These models and their study have been important advances in the understanding of vibration control of construction machinery. Li et al. [16] studied the kinematic accuracy and dynamic performance of spatially deployable mechanisms taking into account the uncertainties of geometric and physical parameters, including the uncertainty of joint clearance. They used the Monte Carlo method to solve the dynamic equation. My et al. [17] proposed a general two-link flexible robot dynamic modeling method, carried out inverse dynamic response analysis, and verified the effectiveness and robustness of the method through numerical calculation. Sun et al. [18] estab-lished a simplified dynamic model of a flexible manipulator with interval gaps and random material characteristics. Their numerical-analysis results showed that joint clearance, component flexibility, and uncertainty all have a strong influence on the kinematic accuracy and dynamic behavior of the flexible manipulator. Meng et al. [19] established a dynamic model of a planar single-link flexible manipulator by using the hypothetical mode method and proposed a fast and stable control strategy based on the system energy to achieve suppression of the flexible rod's vibration. He et al. [20][21][22][23] used the hypothetical mode method, concentrated spring-mass model method, and other methods, including neural network control methods, for flexible link manipulators, to achieve trajectory tracking and vibration suppression. Zhou et al. [24] used space-vector theory and graph theory to study the dynamic coupling effects of multi-arm space robots, focusing on the effects of dynamic parameters, link lengths, and joint variables on coupling. Cuvillon et al. [25] used the theory of modal analysis to derive a dynamic model of the modal space, which is useful for adjusting the damping modal controller; this analysis was verified by experiments on an eight-cable suspension robot prototype.
Researchers have also done considerable work on the chaotic analysis of mechanisms considering the flexibility link, joint clearance, and time-varying stiffness. Chen et al. [26][27][28] established an ideal nonlinear relative dynamics model of the 4-UPS-UPU flexible parallel mechanism and the two-degrees-of-freedom 9-bar mechanism and studied the dynamic response, phase diagram, and Poincaré map of the mechanisms. The mapping and the largest Lyapunov exponent verified the existence of chaos in the complex mechanical system. Cao et al. [29] introduced the fractional-order model of a magnetically coupled broadband energy harvester under low-frequency excitation, using phase trajectories, power spectra, Poincaré maps, and bifurcation diagrams to analyze fractional-order damping, excitation amplitude, and frequency pairs. Lampart et al. [30] discussed the dynamic characteristics of a nonautonomous double pendulum model under soft-stop dual-harmonic excitation, and performed simulations for a specific range of excitation frequencies, confirming that many bifurcations change the characteristics of induced motion into regular, quasi-periodic, and chaotic motions. Jiang et al. [31] studied a cracked gear system in a coal cutter considering the environmental multiple-frequency excitation force, and analyzed the amplitude-frequency relationship of the model's main resonance using the multi-scale method; they also calculated the dynamic response of the model by the incremental harmonic balance method. Qi et al. [32] proposed a spherical chaotic system based on appropriately combining a conservative system with a dissipative system under the celestial mechanism and studied the rich dynamics of four types of spherical system: having spherical four-wing chaotic attractors, spherical ring-like chaotic attractors, spherical periodic orbits, and spherical sink orbits. Based on previous theoretical studies, Moysis et al. [33] proposed an improved HyperJerk system, with a hyperbolic sine function as the only nonlinear term, and studied the path to chaos, antimonotonicity, crisis and coexistence attractors, etc. Miranda-Colorado et al. [34] proposed a hybrid speed controller for flexible joint robots, which can induce chaotic motion on the manipulators of such robots with multiple degrees of freedom. Based on Lyapunov theory, the asymptotic stability of the tracking error signal is verified when the controller is used.
This paper reports our work analyzing the dynamic response and chaotic behavior of a controllable flexible robot. Using the finite element method and Lagrangian equations, an elastodynamic model of the robot is proposed. We establish the law of motion of the robot's end trajectory, compare the natural frequency of the system by using a numerical method and the hammer test method to verify the reliability of the nonlinear dynamic model. The nonlinear dynamic response of the system is analyzed, including the time-domain diagram, the phase diagram, the Poincaré map, the maximum Lyapunov exponent, and the bifurcation diagram. The chaotic state of the system is identified. Finally, the relationship between the motion parameters of the system and the maximum Lyapunov exponent is analyzed, including the radius and angular velocity of the circular trajectory. This lays a foundation for further analysis of the nonlinear characteristics of the system.

Kinematic analysis of the controllable flexible robot
The vibration of an industrial robot will reduce its efficiency, shorten its working life, and reduce its accuracy, so it is of great importance to suppress the robot's vibration. The controllable robot applies the "multidegrees-of-freedom controllable mechanism" theory to the field of robot mechanics. This approach installs the drive motor and reducer on the frame as much as possible and increases the local closed-chain mechanism, thereby reducing the joint's large moment of inertia. The prototype structure consists mainly of a waist, left crank, pull rod, rotating arm, rear arm, auxiliary rod, right crank, forearm, and wrist. The design of the prototype can improve the load capacity and working accuracy while satisfying the requirements of light weight, fast movement, high reliability, and high stability. In kinematic analysis, in order to study the movement of the wrist, one can overlook the rods that make up two parallelograms, and simplify the movement into a classic tandem mechanism, considering only the waist, rear arms, rotating arms, forearms, and wrists. The corresponding joints are all revolving joints, and the angular displacement is used as the relative displacement to meet the requirement that the operation adapts to a complex trajectory. For each joint, the zaxis is in the direction of right-handed rotation, the x-axis is usually in the direction of the common perpendicular to the two z-axis, and the y-axis is always 550 C. Ban et al. Fig. 2 Schematic diagram of a controllable flexible robot perpendicular to both the x-axis and z-axis, as shown in Fig. 2. A corresponding local coordinate system is established at each joint. The rotation axis of joint i is z i−1 , the first coordinate system (coordinate system 0) coincides with the global coordinate system, and the i-th link is located between joints i − 1 and i, in turn analogy. The kinematic relationship between adjacent members of the manipulator can be described by the Denavit-Hartenberg (D-H) [35] model, for which the transformation matrix is where θ i ∈ R n (i = 1, 2, ..., 6) is the generalized coordinate of the joint of the manipulator, d i is the offset of the connecting rod, a i is the length of the connecting rod, and α i is the torsion angle of the connecting rod. The total transformation from the robot-wrist coordinate system to the global coordinate system is 0 Given the angle of the joint θ i , the position and posture of the wrist can be determined by using Eq. (2). The first 3 × 3 columns of the matrix 0 T 6 are the postures of the wrist, and the fourth column represents the position of the wrist. It can be written in the following unified form: where p is the position of the wrist. By differentiating Eq. (3), the relationship between the moving velocityṗ of the wrist and the angular velocity dθ/dt of the joint can be calculated. where For a definite motion ϕ( p), in order to meet the requirements of high-precision control of robot motion, a closed-loop synchronous control scheme is introduced, which includes, for the joint's angular velocity and angular acceleration, where k v and k p are negative definite matrices; e 1 and e 2 are the rigid displacement error and velocity error of the wrist, respectively, andθ s is the initial joint angular velocity.

Elastodynamic model of the controllable flexible robot
Given that modern machinery is developing rapidly in the direction of higher speed, greater precision, and lighter weight, the influence of component deformation on the stability and operating accuracy of the robot cannot be ignored. In order to make the controllable flexible robot achieve better dynamic performance, it is necessary to carry out an elastodynamic analysis. Based on the Lagrangian equation and finite element method, an elastodynamic model of the controllable flexible robot can be established. To simplify the modeling, we make the following assumptions: (1) The waist is treated as a rigid member, the other members as elastic. (2) According to the kinematic assumption, by ignoring the coupling of rigid-body motion and elastic deformation, the real motion can be regarded as a linear superposition of the two. (3) Using the assumption of instantaneous structure at a certain mechanism position, the rotating coordinate system is treated as a static coordinate system, the rotation angle is treated as a constant, and the rotation matrix is treated as a constant matrix. (4) The mass of the bearing mounted at the hinge is treated as a concentrated mass, which is superimposed onto the system mass matrix, ignoring the influence of clearance and friction. Figure 3 shows a typical beam-element structure diagram. The element coordinate system Ax y, which is fixed to the beam element, is established in the diagram. The transverse displacement and longitudinal displacement of the beam element are functions not only of the element coordinate x, but also of time t. It is assumed that the transverse displacement and longitudinal displacement of any point are expressed by w(x, t) and u(x, t), respectively. The transverse displacement function uses a fifth-order polynomial, and the longitudinal displacement uses a first-order polynomial:

Model of beam element
In these equations, δ i is the unit generalized coordinate, and the generalized coordinate array is δ i = [δ 1 , δ 1 , · · · , δ 8 ] T . N i A , N i B are the displacement functions, which depend on the unit coordinate x. and Here e = x/L; it is called the relative coordinate of the beam element, where L is the length of the element.

Kinetic energy of beam element
Assuming that the cross section of the beam element is a uniform rectangle and that the mass of the section is concentrated on the axis, and ignoring the rotational kinetic energy of the section, the kinetic energy T of the beam element can be written where A is the cross-sectional area of the beam element and ρ is the density of the beam.
According to the simplified model assumption (2), it can be concluded that where u r (x, t) and w r (x, t) are the rigid-body displacements of points along the x and y axes, respectively. Bringing the real movement of the mechanism into Eq. (9), the kinetic energy of the beam element can be expressed as where m e is the element quality matrix,

Deformation energy of beam element
When deriving the deformation energy, only the bending deformation of the beam under the action of the bending moment and the tensile and compressive deformation caused by the action of the axial force are taken into account. Other kinds of deformation, such as shear, are ignored. According to the mechanics of materials, the deformation energy can be written where E is the elastic modulus of the beam element and I is the moment of inertia of the cross section of the beam element.
By substituting and simplifying Eq. (13), we get where k e is the element stiffness matrix.

Differential equation of beam element
Incorporating Eqs. (11) and (14) above into the second type of Lagrangian equation, one gets The dynamic equation of the beam element in a rotating coordinate system is where f is the element generalized force array, including the external force array f i and the force array between the units p i , which will cancel each other in the process of mechanism assembly, andδ r is the element rigid-body acceleration array.
In order to transform the assembly of the equation of motion of the elements into the equation of motion of the system, a new generalized coordinate array of elements is introduced: According to the assumption (3), the relationship between δ i and δ ie is where R i is the coordinate transformation matrix. The differential equation of motion of the beam element in the global coordinate system is (20) where m = R T m R , k = R T k R, and f = R T f R are all 8 × 8 symmetric matrices.

Differential equation of the controllable flexible robot
The left crank, right crank, and wrist are cantilever beams. The pull rod and rear arm are each divided into two units. The other members are all single units. A total of 11 units make up the controllable flexible robot system. The generalized coordinate array of the system is defined as where N u = 42 is the total number of generalized coordinates.
The relationship between δ ie and U is written where B i is the coordinate harmony matrix.
In the process of studying the actual working of the controllable mechanism robot, the influence of damping cannot be ignored. In this work, the elastodynamic Table 2 Material characteristic parameters and moment of inertia of the controllable flexible robot    [36][37][38][39]. All of the unit dynamics equations are assembled to obtain the differential equations of motion of the system: whereÜ r is the system rigid-body acceleration array.
F ie are the system mass matrix, damping matrix, stiffness matrix, and generalized force array, respectively, and N e is the total number of elements.

Dynamic response analysis of the controllable flexible robot
The elastodynamic equation of the system, Eq. (23), is a time-varying strongly coupled nonlinear differential equation. In our work, the direct-integration implicit Newmark method and the fourth-order Runge-Kutta algorithm are used to solve the dynamic response of the mechanism. The structure parameters, material characteristic parameters, and moment of inertia of the controllable mechanism robot are shown in Tables 1 and 2 and Table  2. In analyzing the dynamic response of the controllable flexible robot, the period of the motion is defined as T , with T = 2π ω . The integral step is 0.01s. The initial position of the robot wrist is specified as where R is the radius of the circular trajectory.
The controllable mechanism robot moves according to the circular trajectory described by the following equations: Here, ω is the angular velocity of the circular trajectory.
Through numerical simulation, the characteristics of the natural frequency of the system are analyzed, and the time-history diagrams of the elastic displacement, velocity, and acceleration of the robot wrist in the x direction, the y direction, and around the z direction are obtained.

Analysis of frequency characteristics of the controllable flexible robot
The natural frequency of the mechanism is determined by the inherent characteristics of the mechanism (such as mass, shape, and material), and to changes in the position and attitude of the mechanism. By studying the natural frequency of the mechanism, we can grasp the vibration of the system as a whole, avoid a resonance phenomenon of the system, and benefit by the stability of the system. The average value of the natural frequency can express the overall change trend of the mechanism frequency in the entire motion cycle and is an important indicator for evaluating the dynamic characteristics of the mechanism. Among them, the firstorder natural frequency (fundamental frequency) plays an important role in the dynamic response of the system, which is the focus in both engineering applications and theoretical research. By using the classical multi-point excitation singlepoint pick-up modal test method in the force measurement modal test, 14 tested points are selected and marked on the rod mechanism of the controllable flexible robot, and a three-way accelerometer is installed at the end of the robot manipulator, serving as the pick-up point in the modal test. We tap each measuring point in turn, collect the force signal and acceleration signal at each measuring point, then calculate the fundamental frequency of the system. The experimental measurement system for the natural frequency of the controllable flexible robot is shown in Fig. 4, where the robot is stopped in the upright position for testing through the operating handle of the control system. The joint angular position is θ = [0 π/2 0 0 0 0] T . The modal frequency of the system obtained by experiment is 6.836 Hz. The natural frequency, calculated by the numerical simulation method, is 6.84 Hz. The difference between the experimental value and the numerically calculated value is 0.11%, confirming that the elastodynamic model of the controllable flexible robot has high reliability.
Based on the circular motion trajectory given by the mechanism, MATLAB is used to program the numerical simulation. When ω = 0.5π , the relationship between the radii of different circular trajectories and the first-order natural frequency of the system is analyzed, and when R = 0.05m, the relationship between the velocity of different circular trajectories and the first-order natural frequency of the system is analyzed.
It can be seen from Fig. 5 hat the first-order natural frequency of the system decreases as the radius of the circular track increases. With an increase of the angular velocity of the trajectory, the value of the firstorder natural frequency also shows a downward trend, but the local curve fluctuates slightly. Since the natural frequency of each rod is affected by the movement of the mechanism, the overall natural frequency of the mechanism varies with its position.

Vibration response analysis of the controllable flexible robot
According to the law of motion of the circular trajectory of the system, the circular trajectory radius R of the controllable flexible robot is set to 0.05 m, the trajectory angular velocity ω is 0.5π rad/s, and the system material is ordinary carbon steel. The Newmark method is used to solve the nonlinear equations of the system, and the elastic displacement, velocity, acceleration, amplitude spectrum of the robot wrist are recorded. Figures 6, 7, 8 show the elastic displacement, velocity, acceleration, and amplitude spectrum curves of the controllable flexible robot.
As the mechanism follows a circular trajectory, the elastic member of the system will be elastically deformed, causing errors in the vibration of the end manipulator arm, which in turn results in a decrease in the accuracy of the motion of the controllable flexible mechanism. It can be seen that the elastic movement displacement, speed, and acceleration of the mechanism in the x, y directions and around the z direction all have slight oscillations, with no obvious regularity. The frequency domain curves in all directions show the fundamental frequency of the system, which further illustrates the reliability of the elastodynamic modeling of the system.

Chaotic motion analysis of the controllable flexible robot
Although chaotic motion has long-term unpredictability, it does have unique mathematical characteristics that are very important for understanding regular aspects of chaos. In a phase plane diagram with velocity on the ordinate and displacement on the abscissa, one can make a trajectory that changes with time. If a closed trajectory appears in the phase plane, the system has a periodic solution; otherwise, the system may produce chaotic motion. The phase diagram of the controllable flexible robot is shown in Fig. 9. This figure depicts the phase trajectory diagram in the directions of x, y, z, and co-direction. It can be seen from the diagram that the phase trajectory is chaotic and disordered, with the trajectory remaining in the same region.
Through continuous mapping of the system, different forms of phase points or phase trajectories can be observed in the Poincaré map. According to its topological properties, it can be judged whether a system has period 1 motion, period k motion, quasi-periodic motion, or chaotic motion. The isolated points, finite (k) isolated points, closed curves, and uncountable point sets distributed in a certain area appearing on the Poincaré section can, respectively, represent the period 1 motion, period k motion, quasi-periodic motion, and chaotic motion. The Poincaré section of the controllable flexible robot is shown in Fig. 10.
It can be seen from this figure that the Poincaré maps of the robot in the x and, y directions and around the z direction all consist of scattered and disorderly dense points, throughout the graphic area. Figures 9 and 10 both show that the controllable flexible robot system is in a state of chaos.
The Lyapunov exponent is a quantity that describes the degree of divergence of adjacent orbits in the phase space. In this work, the small data volume method [40] is used to estimate the maximum Lyapunov exponent. The concrete recipe for calculation is as follows: where q is the number of nonzero values of d j (i); and use the least squares method to draw a regression line, the slope of which is the maximum Lyapunov exponent.
When the Lyapunov exponent is negative, the mapping converges finally to a certain point; when it is zero, the mapping moves periodically; when it is positive, the mapping is in a chaotic state. The Lyapunov exponent of the controllable flexible robot in the directions of x, y, z and in combination is shown in Fig. 11.
As shown in this figure, the largest Lyapunov exponents are 0.4036 in the x direction, 0.1910 in the y direction, 0.1755 around the z direction, and 0.3058 in the combined direction. The fact that the maximum Lyapunov exponent of the system is greater than zero in all directions confirms that the controllable flexible robot is in a chaotic state. Figure 12 depicts bifurcation diagrams for excitation frequency of the robot in the x direction, the y direction, around the z axis, and in the co-direction. According to this figure, when the external excitation changes from 0 to 0.16 Hz, the motion state of the wrist in all directions changes from periodic motion to chaotic motion. This shows that there is indeed chaos in the elastodynamic response of the controllable flexible robot. In order to study the influence of the circular trajectory radius R on the chaotic motion of the system, the relationship between this radius and the maximum Lyapunov exponent of movement in each direction is established, as shown in Fig. 13.

(d)
As can be seen from this figure, with an increase of the radius R of the circular trajectory, the maximum Lyapunov exponents in the y direction and around the z direction show a downward trend, the maximum Lyapunov exponent in the x direction oscillates around 0.4, while the maximum Lyapunov exponent in the codirection increases at first and then decreases.
6.2 The influence of the angular velocity of the circular trajectory on the chaotic motion of the system In order to study the influence of the angular velocity of the circular trajectory on the chaotic motion of the system, the relationship between the angular velocity of the circular trajectory and the maximum Lyapunov exponent of the movement in each direction was established, as shown in Fig. 14.
It can be seen in this figure that with an increase of the angular velocity ω of the circular trajectory, the maximum Lyapunov exponents in the x and y directions show an overall upward trend. When the angular velocity ω of the circular trajectory is 1 rad/s, the maximum Lyapunov exponent in the combined direction is the largest, 0.6323.

Conclusion
This paper establishes a nonlinear elastic dynamic model of a controllable flexible robot, obtains the nonlinear dynamic response characteristics of this new type of robot system, identifies the chaotic phenomenon in the system, and reveals the reason for the system's abnormal vibration. The main conclusions are the following: (I) Based on the Lagrange equation and the finite element method, the elastodynamic model of the controllable flexible robot system was established. The natural frequency of the robot was tested by a hammering test, and the validity of the nonlinear dynamic model was verified by comparing it with the numerical simulation results. (II) The elastic displacement, velocity, and acceleration of the wrist in accordance with the circular trajectory were analyzed. The simulation results show that the motion of the system has no obvious regularity. By analyzing the system phase diagram, Poincaré map, Lyapunov exponent, and bifurcation diagram, it was further verified that there is chaos in the elastic vibration of the controllable flexible robot under the excitation of inertial force. (III) On this basis, the relationship between the changes in the main motion parameters and the maximum Lyapunov exponent was analyzed comprehensively and systematically. It was revealed that the radius of the trajectory and its angular velocity both have some influence on the nonlinear motion of the controllable flexible robot system. In future work, vibration suppression and joint clearance in controllable flexible robotic systems may be investigated, and this work will also be extended to 3D models.