Delay-induced bifurcations in collocated position control of an elastic arm

The interference of the elasticity of a single robotic arm and the unavoidable time delay of its position control is analysed from nonlinear vibrations viewpoint. The simplified mechanical model of two blocks and a connecting spring considers the first vibration mode of the arm, while the collocated proportional-derivative (PD) control uses the state of the first block only and actuates also there. It is assumed that the relevant nonlinearity is the saturation of the delayed control force. The linear stability analysis proves that stabilizable and non-stabilizable parameter regions follow each other periodically even for large spring stiffnesses and for tiny time delays. Hopf bifurcation calculation is carried out after an infinite-dimensional centre manifold reduction, and closed-form algebraic expressions are given for the amplitudes of the emerging oscillations. These results support the experimental tuning of the control gains since the parameters of the arising and often unexpected self-excited vibrations can serve as a guide for this practical procedure.


Introduction
The position and trajectory control of mechanical systems have been investigated extensively in the literature [27]. There are two basic approaches: the feed forward and the feedback methods. In the feed forward case, the effect of a command is estimated by means of a mechanical model, and the actuation is chosen in a way that the dynamic system performs the desired motion. Input shaping is a widely used feed forward control strategy that was first published in the 50s by Smith [29], and is still a research topic of control strategies [2,6,19,20]. In contrast, the feedback control approach measures the actual states of the system and acts based on the measured states. However, using a feedback loop can be time-and resource-consuming; therefore, the feed forward and feedback approaches are occasionally combined [3].
The control strategies are often designed without the inclusion of time delay [4,14,25,34]. However, due to the finite time of data processing, data transmission and the delayed response of the actuators, a certain amount of time delay always occurs, which may influence the stability properties. This is even more relevant in case of human operated machines with human sensing and actuation [22].
Mathematical modelling of delayed control of mechanical systems has two essential ways. In case of digital control, the states of the system are measured with a given sampling frequency, and the corresponding mathematical model leads to a discrete time system [9,13]. However, when the sampling time is small compared to the data processing and transmission time, continuous time delay description can be applied, which leads to delay differential equations (DDE) [21,32]. Similar DDE models are used in case of human-controlled machines [12].
In many robotic applications, the effectors (like cameras, electrode holders of welding machines or polishers) are located at the end of an arm, the position of which is prescribed by the given technological process. If the corresponding robotic arm is not perfectly rigid, the end effector may oscillate before it reaches its desired position in space since the arm is actuated at its other end as shown in Fig. 1.
There are several solutions for this vibration problem. It is possible to use feedback control in a way that the applied actuator force is determined by the actual state of the end effector. This is called non-collocated control since the sensor of the state of the end effector is at a different location where the actuator is. There are a couple of practical problems with this solution: apart from the increased costs of sensing the position and velocity of the end effector in space, the non-collocated control may be sensitive for stability issues even if the time delays are negligible in the control loop.
The current practical solution is to use collocated position control, where the state of the arm at the actuator is sensed and used in the feedback loop. This way, the position of the end effector is only estimated but the stability of the system can easily be guaranteed for delay-free systems. However, if the delay of the actuation is not negligible, the collocated position control may also become unstable even for quite rigid arms and relatively small time delays. The present study discusses the corresponding conditions of stability and the possibly arising nonlinear vibrations.
Although the delayed PD control and the destabilizing effect of the so-called unmodelled high-frequency dynamics is well-known in the theory of robot control [3], the resolution of the corresponding stability issues and vibration problems is usually left for the experimental gain tuning procedure in practice. An essential nonlinear vibration issue is whether the bifurcations at the stability boundaries are super-or subcritical. In case of delayed control, these boundaries are often subcritical [15,16,18,33], while the saturation of the control force tends to cause stable vibrations around unstable equilibria [8]. The exploration of these phe- In the present case, the delay-induced self-excited vibrations are investigated in collocated position control of an elastic arm with an aim to support, for example, the design of controlled worktables of machine tools carrying elastic workpieces [26] or elastic robotic arms subjected to heavy end effectors [17].
The outline of this paper is as follows. In Sect. 2, the simplified mechanical model is constructed in the presence of collocated feedback control with delay and saturation. The stability of the linearized system is examined and stability charts are constructed in Sect. 3. This is followed by the Hopf bifurcation calculation to characterize the emerging self-excited vibrations at loss of stability. The conclusions are summarized in Sect. 5.

Mechanical model
The mechanical model presented in Fig. 2 is a low degree of freedom approximation of the collocated control of the elastic arm in Fig. 1. Two blocks slide along a smooth horizontal straight path, which are connected with a linear spring of stiffness k that represents the elasticity of the arm in Fig. 1. The displacements of the blocks of masses m 1 and m 2 are x 1 and x 2 , respectively. The generalized coordinates are chosen in a way that the spring is relaxed when x 1 = x 2 . Between the block m 1 and the ground, dry friction is considered, which models the friction that always occurs in the drive chain of the actuator.
In accordance with the collocated control strategy, the position and velocity of mass m 1 are measured and fed back with a proportional-derivative controller. The Denoting the time derivative by dot, the governing equation of the above described system assumes the form where is the simplest model of Coulomb friction force, the maximal absolute value of which is C.
Considering also the saturation as a nonlinearity, the expression of the delayed control force takes the form with the proportional and derivative control gains K p and K d , respectively. In the close neighbourhood of the desired zero position, the control force is assumed to be linear, which leads to an upper estimate of the static position error This error occurs when the system stops at δ], and consequently, the absolute value of the control force |F| becomes less than the Coulomb friction C and the block m 1 sticks to the ground there. This means that the final position error δ is inversely proportional to the gain K p after the control task is completed. Accordingly, the proportional gain should be maximized while keeping the system stable. This is challenging because the presence of time delay always leads to dynamic loss of stability as the proportional gain is increased.
Since the dry friction has a negative power, it has a stabilizing effect on the system [1]. Consequently, when the dry friction is neglected in the following stability and bifurcation analysis, a conservative estimate is applied from engineering viewpoint, while the maximization of the proportional gain still remains an essential goal.
To reduce the number of parameters, let us introduce the dimensionless timet = t/τ and the new dimensionless parameters: which are the mass ratio μ, the dimensionless proportional gain k p , the dimensionless differential gain k d . The parameter α can be considered either as a dimensionless delay or as a dimensionless natural frequency of the uncontrolled mechanical system when the block of mass m 1 is fixed. We also introduce a new saturation parameter Let us calculate the Taylor series of the governing equation up to third order around the equilibrium position x 1 = x 2 = 0, where both blocks are at rest, the spring is relaxed, and the control force is zero. Dropping the tilde, the nonlinear DDE assumes the form: with y = col[x 1ẋ1 x 2ẋ2 ], dot referring to derivative with respect to the dimensionless time, and where

Linear stability analysis
First, let us examine the linear part of the governing equations. With respect to the dimensionless characteristic exponent λ, the corresponding characteristic function yields the quasi-polynomial characteristic equation: which has got infinitely many roots. These characteristic roots determine the stability of the system: it is exponentially stable if and only if all the roots have negative real parts. For the construction of the stability chart, the characteristic roots are checked whether a root crosses the imaginary axis through the origin or a pair of complex conjugate roots enters the right-hand side of the complex plane.
The substitution of λ = 0 into the characteristic equation yields the static section of the D-curves (called static D-curve), where saddle-node bifurcation may occur; this is at k p = 0, along the k d axis in the (k p , k d ) parameter plane (see Fig. 3).
The so-called dynamic D-curves related to possible Hopf bifurcations can be determined by substituting λ = iω into Eq. (12) where ω = 0 is the dimensionless angular frequency of the emerging selfexcited vibrations. Following the steps of the so-called D-subdivision method [11], the separation of the real and imaginary parts of the complex characteristic equation can be rearranged with respect to the critical control gains as a function of the parameter ω:  (13) and (14) are singular at ω = α, and for α < π/2, the left-hand side and right-hand side limits are plus and minus infinity, respectively. The dynamic D-curve shown in Fig. 3 starts from the origin of the (k p , k d ) parameter plane, crosses the origin again at ω = α √ 1 + μ, and spirals outwards counter clockwise.
The next important property of the system is the so-called root tendency, that is, in which direction the characteristic roots cross the imaginary axis as parameters vary. This can be determined with the help of implicit differentiation of the characteristic Eq. (12) with respect to the gain k p where the characteristic exponent is considered as a function λ(k p ) of this parameter: where prime denotes the derivative with respect to k p , and The denominator of Eq. (15) is positive for ω > 0 and ω = α. Thus, the numerator of (15) specifies the sign, and so the direction of the root tendency as the dimensionless control gain k p increases. With the help of this, the number of the unstable characteristic roots can be determined in the disjunct parameter domains of the stability chart (see Fig. 3), which are defined by the stability boundaries. Figure 3 presents a D-shaped stable region. The system may loss its stability in a static way through the saddle-node bifurcation or in a dynamic way through the Hopf bifurcation where the emerging self-excited vibration has the dimensionless angular frequency ω ∈ [α √ 1 + μ, π/2]. The stable area shrinks as the dimensionless delay α increases, and it disappears at: it is a critical value for α, that is, for the time delay of the system. However, it is not an absolute maximal delay where the system is still stabilizable. If we increase α further, there will be further intervals for α where the system can be stabilized again (see Fig. 4) with appropriate control parameters. In this context, the linear system with fixed parameters will be called not stabilizable if no proportional and derivative control gains exist to stabilize it. The existence of the stable region is related to the slope of the dynamic D-curve at the origin. The initial slope of this boundary at ω = 0 is independent of the dimensionless delay α: However, when this boundary curve crosses the origin again at ω = α √ 1 + μ, the same derivative is which already depends on α. Clearly, a critical value of α occurs when this derivative changes sign, while it tends to ±∞. The first such critical value is at α cr,1 as given in (18). As the delay increases further, that is, α > α cr,1 , the first branch of the dynamic D-curve for ω ∈ [0, α) bends to the left, and in the same time, the second branch for ω ∈ (α, ∞) has a tangent line at ω = α √ 1 + μ that rotates counter clockwise according to Eq. (20). Thus, the second branch of the dynamic D-curve intersects the D-shaped region defined by the first branch for ω ∈ [0, π/2] as a "windscreen wiper" turning it again and again to stable and unsta- It is interesting to compare these stability charts to the ones of the PD-controlled single degree of freedom (DoF) system, where k → ∞ and m = m 1 + m 2 . The corresponding dynamic D-curves are presented with the red dashed lines in Fig. 4. As α increases, the first branch of the dynamic D-curve of the original system tends to this Hopf-type stability boundary of the single DoF system, which yields that the reappearing Dshaped stable region approximates the stable domain of the single DoF system (see for example the panel of Fig. 4 with α = 6). This agrees with our intuition, since, by definition, α 2 is directly proportional to the stiffness of the spring, thus, the connection between the blocks becomes rigid as α → ∞. Note that while these curves practically coincide for α > 6, they are not necessarily stability boundaries of the two DoF system due to the above described periodic "windscreen wiper effect" caused by the rotating tail of the second branch of the dynamic D-curve for ω ∈ (α, ∞).
The lower panel of Fig. 5 presents the stability boundaries in the dimensionless parameter plane (α, k d ) at fixed values of the dimensionless control gain k p . One can observe that increasing α from 0, the stable domain shrinks, disappears and when it appears again, it is even larger than in the case of α = 0. Note, however, that this phenomenon occurs only in the plane of the dimensionless control parameters. If one transforms them back to K p and K d , then at α = 0, that is at τ = 0, the whole positive quadrant of the (K p , K d ) parameter plane is stable, while above α cr,2 , the reappearing stable domain is bounded and its size is smaller.
A practically relevant consequence of the dynamic analysis is related to the calculation of the possible maximal proportional control gains K p,max that ensures the smallest possible static position error δ min = C/K p,max in (4). Panel (a) of Fig. 6 presents the maximal dimensionless proportional gains k p,max against the dimensionless delay α, which are determined by means of the rightmost parameter points of the stable domains (if exist) in the stability charts of Fig. 4. It can be observed that k p,max is even larger in the second domain of stability than in the first. However, the dimensional proportional gain K p,max = k p,max m 1 /τ 2 may have very large values in the first stable domain  Fig. 6) because τ → 0 yields K p,max → ∞.
In the meantime, since the dimensionless parameter α also includes the spring stiffness k (see (5)), one can calculate K p,max also as a function of the spring stiffness [see the parameter point P in panel (c) of Fig. 6]. This results in the overall smallest static position error in the second stable domain at It clearly shows that decreasing the time delay improves the positional accuracy substantially.

Hopf bifurcation calculation
As a parameter point leaves the stable region of the (k p , k d ) plane at the stability boundary presented by the dynamic D-curves, a limit cycle may emerge, the amplitude of which can be approximated with the help of the Hopf bifurcation analysis. Its first step is to carry out the so-called centre manifold reduction, that is, to project the system from the infinite-dimensional state space of the DDE onto the two-dimensional invariant subspace that is tangent to the plane spanned by the critical eigenvectors [11,24]. Introduce the function x t : R → X R 4 with the shift x t (ϑ) = x(t + ϑ), with ϑ ∈ [−1, 0], where the minimum of ϑ corresponds to the length of the delay, which is 1 in the dimensionless form, and X R 4 is the space of continuous functions mapping the interval [−1, 0] into R 4 . Furthermore, one can formulate the delay-differential equation (7) as an operatordifferential equation (OpDE) [10]: where the operators are defined as Here, L and R are the coefficient matrices of the nondelayed and delayed state vectors in Eq. (7), respectively. The characteristic exponents of the linear part of (7) are the same as those of the linear operator A [11].
For the centre manifold reduction, we need the right eigenvectors of A corresponding to the critical eigenvalues λ cr = ±iω. Then, the real eigenvectors s 1,2 ∈ X R 4 satisfy (see [30]): which is a boundary value problem for s 1,2 (ϑ) according to the definition of the linear operator A in Eq. (23). The solution assumes the form s 1 (ϑ) = S 1 cos(ωϑ) − S 2 sin(ωϑ), where S 1,2 ∈ R 4 satisfy the linear homogeneous algebraic equation as it follows from the corresponding boundary conditions. Let us choose the first components of S 1 and S 2 to be 1 and 0, respectively; then, (27) yields In order to project the system onto the plane spanned by the eigenvectors s 1 and s 2 , one needs the left eigenvectors n 1,2 of the linear operator A associated with the critical eigenvalues ∓iω.
Introduce the adjoint operator of A by where denotes the transposed conjugate matrix and also the adjoint operator. The real eigenvectors n 1,2 ∈ X R 4 satisfy From which with N 1,2 ∈ R 4 satisfying L + R cos ω − (ωI + R sin ω) ωI + R sin ω L + R cos ω The solutions include the two free parameters N 12 and N 22 that are determined by means of the orthonormality conditions n 1 , s 1 = 1, n 1 , s 2 = 0, where ., . denotes the inner product. With the definition of the inner product (see "Appendix A"), (35) can be expanded as which yields Here, a and b are expressed in (16) and (17), respectively, which also appear in the root tendency Eq. (15). In order to restrict the dynamics to the invariant centre manifold, the infinite-dimensional state x t is projected to the right eigenvectors s 1 and s 2 by means of the corresponding inner product between the left eigenvectors n 1,2 and x t . The infinite-dimensional part w of the state x t is calculated by simple subtraction. Accordingly, the new scalar coordinates z 1,2 on the centre manifold are introduced by and With these new variables, the operator differential Eq. (22) can be reformulated as ⎡ where o : R → X R 4 is a zero operator, O : X R 4 → R is a zero functional; for details, see "Appendix B".
In order to get the first Fourier term of the oscillation in the centre manifold, the first two rows of Eq. (42) should be calculated up to third order: where the indices k and l are nonnegative integers. Due to the symmetry of the nonlinearity, the nonlinear operator F does not have second-order terms, which implies that w is at least O(y 3 1,2 ), and so it effects only the fifthand higher-order terms in (43).
The so-called Poincaré-Lyapunov coefficient Δ can be calculated with the help of the Bautin formula (see "Appendix C"). Since one obtains f 12 = g 21 = g 03 = 0, This is always negative, that is, the Hopf bifurcation is always supercritical along the dynamic D-curves of Fig. 4.
In the centre manifold, the oscillation close to the Hopf bifurcation boundary can be approximated as where the amplitude A can be expressed as This yields that for proportional gain k p close enough to its critical value k p,cr , the self-excited vibration can be approximated as where S 1,2 are given in (28). This way, the substitution of (15) and (44) and These algebraic expressions fit well to the amplitudes created with DDE-BIFTOOL [5,28] for the original nonlinear system given in (7)-(10) (see Fig. 7). Clearly, the higher-order nonlinear terms result in an increasing deviation between the analytical and numerical results for k p > 1.1 k p,cr .

Conclusion
The delayed position control of two blocks connected through a linear spring was investigated as a simplified model of an elastic arm. A collocated proportionalderivative control force was considered with saturation nonlinearity. The linear stability analysis shows that the stable region in the plane of the dimensionless control parameters (k p , k d ) shrinks and disappears as the time delay increases. Further increase in the delay results in the reappearance of the stable parameter domain, which then disappears and reappears periodically. Note that these stable regions become negligible for the original control gains K p , K d as it follows from the definitions of the dimensionless parameters in (5) for large delays.
However, the interpretation of the same dimensionless stability charts is different from the viewpoint of the stiffness of the linear spring. The increase in this stiffness also increases the dimensionless parameter α [see (5)], which means that in case of a fixed time delay τ , the phenomenon of the periodic appearance/disappearance of the stable control parameter regions is also present as the stiffness parameter varies. These stable domains have about the same size for the original control parameters K p , K d when the delay τ is fixed, and they get close to the D-shaped stable domain of the PD-controlled single DoF system, where k → ∞ and m = m 1 + m 2 .
Also in case of fixed delay τ , the static positional accuracy can be optimized with the help of (21) by tuning the system and control parameters accordingly.
If the stability chart of only the single DoF model is used in the tuning procedure of the control parameters, unexpected instabilities may occur even for relatively large spring stiffness values [see also panel (c) of Fig. 6] where a corresponding robotic arm might be assumed to be practically rigid (see Fig. 1).
The Hopf bifurcation calculation was executed after an infinite-dimensional centre manifold reduction. The Hopf bifurcation is supercritical all along the dynamic stability boundaries. The closed form algebraic expressions for the amplitudes of the arising self-excited oscillations are also checked by the freely available MAT-LAB software DDE-BIFTOOL. The comparison of the calculated and measured vibration amplitudes may help to estimate how far the stability boundaries are, in which direction the control parameter tuning procedure should be continued.
At the boundary of the stable area, the dynamic D-curve can intersect itself, which may correspond to codimension-2 Hopf bifurcation, where two pairs of complex-conjugate characteristic exponents cross the imaginary axis. At this parameter setting, the selfexcited vibrations of the system may become quasiperiodic oscillations containing two distinct frequencies. See, for example, the lower left panel of Fig. 4 with α = 3.5; here, the Hopf-Hopf point is denoted with H 2 , and the two independent dimensionless angular frequencies are ω 1 = 1.029 and ω 2 = 4.578. The codimension-2 Hopf bifurcations can be analysed with the algorithm presented in [23,31]. The existence of quasiperiodic oscillations is expected from the topological structure of the H 2 bifurcations [7], while its rigorous proof needs further work similarly to the analysis of higher DoF models. Note that the occurrence of quasiperiodic oscillations during the experimental tuning of the control parameters may provide an essential validation of the time-delayed mechanical model used in this study.
The methodology presented here can be generalized for multi-degree of freedom robotic arms, but due to the mathematical complexity, no analytical results are expected. In the meantime, the arising nonlinear vibrations will be qualitatively analogous to the ones observed in the low-dimensional case; its justification is part of a future research.
Funding Open access funding provided by Budapest University of Technology and Economics.

Data Availability Statement
Data sharing is not applicable to this article as no datasets were generated or analysed during the current study.

Conflicts of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.