Bistability and delayed acceleration feedback control analytical study of collocated and non-collocated cases

Stability and bifurcation analysis of a non-rigid robotic arm controlled with a time-delayed acceleration feedback loop is addressed in this work. The study aims at revealing the dynamical mechanisms leading to the appearance of limit cycle oscillations existing in the stable region of the trivial solution of the system, which is related to the combined dynamics of the robot control and its structural nonlinearities. An analytical study of the bifurcations occurring at the loss of stability illustrates that, in general, hardening structural nonlinearities at the joint promote a subcritical character of the bifurcations. Consequently, limit cycle oscillations are generated within the stable region of the trivial solution. A nonlinear control force is then developed to enforce the supercriticality of the bifurcations. Results illustrate that this strategy enables to partially eliminate limit cycle oscillations coexisting with the stable trivial solution. The mechanical system is analysed in a collocated and a non-collocated configuration, depending on the position of the sensor.


Introduction
Robots were originally developed for pick-and-place operations [25], for which extreme geometrical accuracy is not a particular need. Later, industrial robots were also used as a milling platform, first for prototype manufacturing of clay [38,39] mostly for industrial design cases. The repeatability of positioning is more important for machining than for pick-and-place operations, which was enough for the mentioned clay prototype sculpture manufacturing. By the time, geometrical accuracy has increased by a factor of 10 from mm to tenth of mm by applying better sensors and proportional-derivative feedback control. Prototype manufacturing is slowly taken over by plastic 3D rapid prototyping [26,44]; however, for some cases, clay still serves as cheap and durable solution to test design issues of e.g. sculpture-like car chassis prototypes. Indeed, the most common open kinematical chain robots are excellent to manufacture sculpturelike surfaces [9][10][11] due to their workspace, but this versatile configuration lacks rigidity.
Machine tools are created to perform machining on metal parts; their superstructure is designed to serve stiffness and rigidity [29], required to provide the necessary geometrical accuracy. This leads to a design very different with respect to robots, resulting in heavy, large, sturdy and consequently expensive superstructures for machine tools.
Machining of metal parts usually requires high shape accuracy; at first, the reflected stiffness must overcome the corresponding so-called cutting stiffness [43] ensuring the minimum requirement for a stable operation. This is still very much a must even for roughing operations, where accuracy is not an important factor, because of the great forces transmitted during cutting. This requirement is especially hard to achieve by industrial robots, which are generally slender and not particularly stiff [31,32,42]. There are stiffer geometric arrangements of open kinematic chains of robotic arms; however, they generally provide a drastically reduced utilisable workspace, if the cutting is possible at all.
Theoretically, an online control solution is possible to upgrade the robot dynamic performance by built-in velocity and acceleration feedback controllers. On the one hand, robot manufacturers are not taking the risk to open their control development kit to that level. On the other hand, possible factory level deployment of upgraded industrial robots can be accepted by owners either using off-line techniques or utilising a certified built-in control.
Since industrial robot users are also hesitant to accept mounted passive and semi-active embedded solutions to increase the dynamic reflected stiffness of the end-effector, another relatively cheap solution could be using simply feedback sensor signals in the built-in position control of stock industrial robots. This might sound impossible, but stable manufacturing might be achieved with careful parameter management, especially with the synchronisation of delays in the built-in and new feedback loop, most probably using acceleration signals.
Usually, a proportional-derivative control scheme is used in the inner controller of the robot [47] to decrease the positioning error. An additional acceleration feedback controller, which is fed through the conventional inner controller of the robot, is sometimes utilised to mitigate the vibrations arising from the machining processes [7,18,34,48,49] or to compensate for the flexibility of the arm [12,19]. Acceleration feedback is particularly appealing for industrial applications thanks to its reasonable price, compared to position and velocity measurement devices required to measure the state of the robot's end-effector [14,15]. Typically, sensors measure only the state of the motor's rotor, while only relatively recently, sensors for directly (or indirectly) measuring the position of the arm have been proposed [3,36]. Because of the flexibility of the robotic arm, especially at the joints [1], the correct positioning of the accelerometer can have a significant effect on the stability and effectiveness of the control scheme [13].
The system is named collocated if the sensors are positioned near the actuator and non-collocated if they are positioned far from the actuator [2,22]. The two cases are illustrated in Fig. 1.
Feedback controllers usually introduce time delays in the system [41], which is another element potentially having a substantial effect on the stability and effectiveness of the controller, even if it is relatively small [17,23]. Nevertheless, its effect is not necessarily detrimental; indeed, several studies exploit time delay to enhance the stability performance of a controller [5]. For instance, in [4] and [45] time delay is purposely included in the system in order to enhance damping because of the phase shift it generates between the acceleration and the control force. A similar study is performed on a non-collocated system in [37], where the damping properties of the delayed acceleration feedback controller are evaluated.
A usual common experience on slender open kinematic chain robots is that they often encounter structural nonlinearities that could induce unstable limit cycles even in case of achieved stable stationary cutting domain. Especially the joints and the drive of the robotic arm can introduce nonlinearities to the stiffness characteristic, which is another feature relevant for the system dynamics often neglected [31,40]. This aspect was partially studied, for instance, in [8], where a lumped mass model with a cubic nonlinearity controlled by a delayed acceleration feedback signal was investigated; the implications of 1:1 resonances and 1/3 subharmonic resonances on stability were analysed. It should be emphasised that the nonlinearity might limit the robustness of the equilibrium position to external perturbations even if the equilibrium is stable [21].
This study was motivated by the industrial experience on appearing structural nonlinearities and the possible solution to use an acceleration feedback control scheme in the built-in robot controller. The main objective of this work is to analyse the influence of the nonlinear behaviour at the joint of the robotic arm on its dynamics. In particular, we aim at unveiling how it affects the bifurcation characteristic at the loss of stability and to which extent this can be controlled.
We attempted to create the simplest possible equivalent model that is able to produce nonlinear behaviour and mimic the robot position control. However, it should be noted that it is far from the real-life behaviour of robotic arms. A two-degree-of-freedom (DoF) lumped model of a linear robotic arm, subject to an acceleration feedback control, is considered. Two configurations are studied, namely a collocated and a non-collocated one, depending on the position of the accelerometer. First, the stability of both configurations is analytically studied, and then, the effect of the various system parameters on the stability is analysed. The stability of the delayed and non-delayed systems is compared, which highlights the importance of considering time delay in such systems. Later, exploiting a multiple-scale perturbation technique [35], the bifurcations occurring at the loss of stability are analytically studied, providing a partial understanding of the global dynamics of the system. Finally, a nonlinear term is included in the control loop to modify the character of the bifurcation in the case of subcritical bifurcations, with the objective of improving the robustness of the stable solution. Limitations and advantages of the various configurations and control laws considered are carefully analysed.

Mechanical model
The adopted two-DoF mechanical system is shown in Fig. 2.
The model consists of two lumped masses m 1 and m 2 connected by a nonlinear spring k and a linear damping c. The prescribed reference trajectory u d is programmed, such that in ideal circumstances an identical constrained motion u r is imposed to the mass m 1 via a linear spring k 1 with time delay τ r . This gives the dynamical system governed by the following equations of motion where u 1 and u 2 are the displacements of the lumped masses. Prime ( ) indicates derivation with respect to time t. Furthermore, Δu := u 2 − u 1 and the nonlinear spring is considered with a second-order nonlinearity in the stiffness characteristic k(Δu) = k 2 + k nl (Δu) 2 , which will lead to cubic nonlinearities in the system. Apart from the prescribed reference trajectory u d , an additional feedback signal u f , with third-order nonlinear acceleration signal is fed back with time delay τ f and summarised with ideal constrained motion u r . This gives a final constrained motion u , where τ = τ r + τ f is the sum of the delays. The linear and nonlinear acceleration control gains are a 1 , a 2 , a nl,1 and a nl,2 , respectively. The values m 1 , m 2 , c, k 1 and k 2 are assumed to be positive real numbers, while k nl , a 1 , a 2 , a nl,1 and a nl,2 are generic real numbers. Generally, the control force is calculated based on the value of the acceleration of a specific point; therefore, we assume that either a 1 = 0 and a nl,1 = 0 or a 2 = 0 and a nl,2 = 0. These two cases correspond to a noncollocated or collocated controller, respectively. With these considerations, the equations of motion can be written as where . In order to focus on the instabilities induced by the acceleration feedback controller, constant prescribed reference trajectory is assumed With the introduction of perturbations u 1 = u 1 + x 1 and u 2 = u 2 + x 2 , the equations of motion can be written as In order to reduce the number of parameters, we divide both equations by m 1 ; we introduce the dimensionless time tω 1 → t and displacements where ω 2 1 := k 1 /m 1 and r := m 2 /m 1 . We thus reduce the system tö where over-dot (˙) indicates derivation with respect to the new, rescaled dimensionless time t, whileẍ 1τ = x 1 (t − τ ) ,ẍ 2τ =ẍ 2 (t − τ ), with the dimensionless time delay τ ω 1 → τ . Besides, the dimensionless standalone natural frequency ratio γ := ω 2 /ω 1 , the standalone damping ratio χ := c/(2m 2 ω 2 ), with ω 2 2 := k 2 /m 2 , the rescaled control parameters κ 1 := ω 2 1 a 1 , κ 2 := ω 2 1 a 2 , the rescaled nonlinear control parameters κ nl,1 := ω 6 1 a nl,1 rk 1 /|k nl |, κ nl,2 := ω 6 1 a nl,2 rk 1 /|k nl | and μ := sign(k nl ) are introduced.

Stability analysis
In this section, we investigate the linear stability of the system's trivial solution. Firstly, we study the system with a traditional approach, neglecting the time delay. Then, we analyse the delayed system and perform a parametric analysis to evaluate the sensitivity of the stable region to the parameter values.

Non-delayed case
In order to analyse the importance of considering the time delay of the control force, in this section, we study the stability of the collocated and non-collocated systems neglecting the time delay. Imposing τ = 0, regarded to the ideal non-delayed control, the system in Eq. (4) is reduced tö This means that the control actions affect the inertia, while κ nl,1 and κ nl,2 introduces conventional nonlinear terms. By solving the system with respect toẍ 1 and x 2 and linearising it around its trivial equilibrium, we obtain the linear systeṁ where The stability of this system can be investigated with the Routh-Hurwitz criterion [46]. For this purpose, we calculate the characteristic equation which can be rewritten as (9) where I is the identity matrix. Considering the following Hurwitz matrix the system is stable if all the leading principal minors of H are positive. Conditions providing stability of the trivial solution for the collocated and non-collocated cases are discussed below.

Collocated case
If κ 2 = 0, then the four conditions from the Hurwitz matrix are These conditions can be reduced to two sets of conditions The frequency ratio γ and mass ratio r have to be positive real values. Therefore, the following conditions must be fulfilled for the trivial solution to be stable: This result seems in disagreement with respect to what obtained regarding the stability analysis for the delayed case, illustrated in the following section ( Fig. 3(i)). This discrepancy is discussed in Sect. 3.2.

Non-collocated case
If κ 1 = 0, then the four conditions from the Hurwitz matrix are These conditions can be reduced to two sets of inequalities: Since the frequency and mass ratio must be positive real numbers, these conditions are reduced to This result is fully confirmed by the stability charts plotted in Fig. 3(ii). In fact, for all investigated cases, as τ → 0 the system is stable if κ 2 ∈ (−∞, r ).

Delayed case
In order to study the stability of the trivial solution of the system, we temporarily neglect nonlinear terms and we consider the trial solution x 1 = ξ 1 e λt and x 2 = ξ 2 e λt , obtaining By applying the D-subdivision method, when λ = iω, the characteristic function has the following form whose solution provides boundaries which separate regions of the parameter space in which the system has different number of eigenvalues having positive real part. Some of these boundaries mark the stable region of the system.
Separating real and imaginary parts of Eq. (18), we obtain Summing up the squares of these two equations leads to Equation (20) can be directly solved with respect to κ 1 or κ 2 . In particular, we first impose κ 2 = 0, which refers to the case of a collocated controller, obtaining Instead, imposing κ 1 = 0 and therefore referring to the case of a non-collocated controller, we have Equations (21) and (22) will be directly used for defining the stability boundary in the τ, κ 1 and τ, κ 2 spaces for the collocated and non-collocated cases, respectively. Solving Eq. (19) with respect to sin τ ω and cos τ ω, we obtain where From Eq. (23), we can directly obtain the value of τ as a function of ω, i.e.
Equations (21), (22) and (25) enable us to plot the stability diagram in the (τ, κ 1 ) and (τ, κ 2 ) spaces. The diagrams, illustrated in Fig. 3, are obtained spanning a large range of ω values adopting various n values in Eq. (25). The effect of parameter values on stability is discussed in Sect. 3.3. Figure 3(i) suggests that a necessary condition to have stability is that −1 < κ 1 < 1, for the collocated case. This condition can be obtained analytically considering the limit of κ 1 for ω → ∞ in Eq. (21), which gives ±1. This result is typical for neutral delay differential equations as discussed, for instance, in [30]. Interestingly, this conclusion implies a discontinuity with respect to the non-delayed case for τ ≈ 0 and κ 1 < −1. In fact, as demonstrated in Sect. 3.1, the trivial solution of the collocated system is stable for χ > 0 and κ 1 < 1, therefore, also for κ 1 < −1. Figure 3(i) clearly shows that for τ → 0 the system is unstable if Let us investigate the eigenvalues of the system in correspondence with the discontinuity. Figure 4a, b depicts some of the eigenvalues of the system for a very small time delay (τ = 10 −5 ), for κ 1 = −0.9999 and κ 1 = −1.0001, respectively. The four eigenvalues in the insets of the two plots exist for both the delayed and the non-delayed cases. In both cases, they remain in the left part of the complex plane, which suggests that the instability is not related to their perturbations. The two figures show other eigenvalues having the absolute value of their real part approximately 10. However, while for κ 1 = −0.9999 (Fig. 4a) the real part is negative, for κ 1 = −1.0001 (Fig. 4b) it is positive, which confirms the stable and unstable behaviour of the two cases, respectively. Nevertheless, since these eigenvalues are far apart despite having very similar parameter values, they confirm the discontinuity of the spectrum also for κ 1 = −1 and τ = 0. We speculate that the eigenvalues might move from the left to the right side of the complex plane through ±∞ of the imaginary axis, which is consistent with the limit of κ 1 for ω → ∞ in Eq. (21), and with the large value of their imaginary part. Indeed, a similar scenario might be responsible for the spectral discontinuity occurring at τ = 0 and κ 1 < −1.
Previous studies [24,27,28,33,50] illustrate that similar discontinuities in the spectrum are rather common for systems governed by neutral delayed differential equations. For instance, in [28,50], analogous discontinuities are encountered in the delayed control of an elastic beam, both for τ = 0 and for τ = 0. The possible existence of such discontinuities was mathematically proven in [24]. Additionally, Insperger et al. [27] encountered a similar discontinuity for advanced delay differential equations as well. Figure 4c, d provides time series of the system dynamics in correspondence of the discontinuity. Figure 4c shows that, for the non-delayed case, the system is stable and has practically unvaried behaviour for κ 1 slightly larger or smaller than −1. On the contrary, Fig. 4d illustrates that, for the delayed case (τ = 10 −5 ), Fig. 4 (a, b) Eigenvalues of the collocated system for r = 1, γ = 1, χ = 0.1, τ = 10 −5 , κ 1 = −0.9999 (a) and κ 1 = −1.0001 (b); the insets show four eigenvalues of the spectra, which coincide with the eigenvalues of the non-delayed system. c, d Numerical simulations of the system dynamics for the same parameter values; c non-delayed case (τ = 0), d delayed case (τ = 10 −5 ); κ 1 as indicated in the figure if κ 1 < −1, the system rapidly experiences instabilities that cause high-frequency oscillations of increasing amplitude.

Parameter analysis
In this section, we investigate the sensitivity and behaviour of the stable region with changing parameters. The stable region is studied in the plane of the time delay τ and control parameters κ 1 , κ 2 . Generally, one would think that lower time delays would be beneficial for stability. However, the control of the robot and the acceleration feedback controller can introduce relatively large time delays in the system, so the stable region for larger τ values has to be studied. Although the actuator's saturation is neglected in this study, it should be noted that technological constraints limit the value of the control parameters κ 1 and κ 2 ; therefore, they cannot have arbitrarily large values. The parameter analysis of the system is illustrated in Fig. 3.
The light green line in Fig. 3 shows the controller with the best performance, i.e. where the largest characteristic root is minimum for one particular parameter set. It should be noted that this is very close to the boundary of stability; therefore, when one wants to operate the controller close to this line, the effect of the nonlinear terms can be particularly relevant for the system dynamics.

Variations of the mass ratio r
Stability charts for the collocated and non-collocated cases for various mass ratio r values, with the other parameters kept constant γ = 1, χ = 0.1, are displayed in Fig. 3a. Regarding the collocated case (Fig. 3a, (i), the stable region is not very sensitive to variations of the mass ratio r . The positive stable region shrinks with increasing τ values. On the contrary, the negative stable region undergoes only a slight contraction for all the considered τ range. This feature, for the collocated case, will also be confirmed for other parameter values investigated below.
Conversely, the non-collocated system is much more sensitive to variations of r . In general, the stable region enlarges with r , particularly for small τ values. This effect is valid for both positive and negative ranges of κ 2 . We notice that for τ → 0, the upper stability boundary tends to r , as calculated for the non-delayed system in Sect. 3.1.2. Therefore, this parameter is important because it defines the largest positive stable value of κ 2 . In the negative region, as τ tends to zero, the stable values of κ 2 tend to negative infinity. With increasing time delay, the positive stable region shrinks much faster than the stable negative region, and it gets very small with large delays, for all considered r parameters. Interestingly, in the non-collocated case, the stable region is larger for κ 2 < 0 than for κ 2 > 0, being even unbounded for small τ values. On the contrary, for the collocated case, the stable region is generally smaller and bounded for any τ value.

Variations of the standalone frequency ratio γ
Stability charts for the collocated and non-collocated cases for various γ values, with the other parameters kept constant r = 1, χ = 0.1, are represented in Fig. 3b. For both collocated and non-collocated cases, the stable region is larger for smaller γ values. However, in the non-collocated case, the system is much more sensitive to variations of γ . Regarding the collocated case, with larger γ values, the positive stable region abruptly shrinks for relatively small τ values. Considering the non-collocated case, we notice that for relatively small γ values (in Fig. 3b (ii), γ = 0.5 and γ = 0.75), a sort of peak can be given for the maximum positive value of κ 2 . Smaller γ values mean that the system is tuned for larger ω 1 values than ω 2 values, so the robot is stiffer close to the actuator. As expected, results illustrate that the stable region is larger in this case.

Variations of the standalone damping factor χ
Stability charts for the collocated and non-collocated cases for various standalone damping factor χ values, with the other parameters kept constant r = 1, γ = 1, are displayed in Fig. 3c. As expected, the stable region becomes smaller for decreasing χ values in both cases. The shrinkage is more pronounced for larger values of time delay τ . Considering the collocated case, we notice that for sufficiently large χ values, the destabilising effect of the time delay becomes negligible, and the maximum and minimum values of κ 1 which guarantee stability are constant with respect to τ for a large range (see, for instance, the blue curve in Fig. 3c (i)). Furthermore, the negative region does not seem to be particularly affected by small values of χ ; for instance, reducing χ from 0.1 to 0.01 practically leaves the negative stable region of the collocated case unchanged. Regarding the non-collocated case, the destabilising effect of small values of χ is relevant for both positive and negative half-plane of the stability diagram. However, the positive critical value of κ 2 for τ → 0 is unaffected by variations of χ , which means that in the non-collocated case the stable region is much larger than in the collocated case as τ → 0. Nevertheless, even relatively small nonzero τ values shrink the stable region more significantly in the non-collocated case than in the collocated one.
We remark that, although the stable region significantly shrinks in some cases, the system is always stable for zero control force. This means that we cannot talk about a "more stable" configuration. Conversely, we can state that the system is more sensitive to large control gains in some cases. For instance, let us consider the case of τ → 0 and large mass ratio (m 2 m 1 ), which makes the non-collocated configuration stable for very large κ 2 values. This result does not mean that the system is very stable but that a control force applied on the first mass-which is relatively small-proportional to the acceleration of the second mass, has a negligible effect on the system stability, which indeed sounds physically reasonable.

Bifurcation analysis
We now aim at analysing the Andronov-Hopf bifurcations, identified by the D-curves, occurring at the loss of stability, and the effect of the sign of the system nonlinearity on it. In particular, a subcritical bifurcation might limit the robustness of the stable trivial solution in the proximity of the stability boundary, which can cause the onset of self-excited vibrations within the stable region. The method of multiple scales for timedelayed systems is adopted for such analysis [20,35].
First, we transform the system in Eq. (4) into firstorder form, i.e. y = Ly − κ 1 R 1ẏτ + f 1 (26) in the case of a collocated controller, oṙ for a non-collocated one, where The adopted bifurcation parameters are the rescaled control parameters κ 1 and κ 2 .
Next, in order to eventually obtain a convergent solution, we expand the solution of Eq. (26) (or Eq. (27)) as where T 0 = t, T 1 = εt and ε is a small bookkeeping parameter. The derivative with respect to t is trans- while delayed terms are expressed as which upon expansion for small ε becomes The bifurcation parameters κ 1 and κ 2 are rewritten as κ 1 = κ 1,cr + εδ 1 and κ 2 = κ 2,cr + εδ 2 , respectively, where κ 1,cr and κ 2,cr correspond to the κ 1 and κ 2 values at the loss of stability, while δ 1 indicates the closeness of κ 1 to κ 1,cr , δ 2 has an analogous meaning. As explained later, this expansion enables us to have a non-hyperbolic fixed point for order ε 1/2 and secular terms generated for order ε 3/2 . Indeed, this choice of the order of ε is suitable for systems having only cubic nonlinearities, like the one under study. For a more comprehensive explanation of the formalism adopted, we address the interested reader to [35]. Substituting Eqs. (29)-(32) into Eq. (26), we obtain where h.o.t. indicates terms with values of ε higher than order 3/2, y 1τ = y 1 (T 0 − τ, T 1 ), y 2τ = y 2 (T 0 − τ, T 1 ), y 11 , y 21 , y 31 and y 41 are the elements of vector y 1 . In the non-collocated case, κ 1,cr , δ 1 , R 1 andf 1 should be substituted with κ 2,cr , δ 2 , R 2 andf 2 . For the sake of brevity, in the following, equations related to the noncollocated case will be omitted unless necessary. Separating terms of different order of ε, we obtain Order ε 1/2 Equation (34) consists of the linearised system of equation for a critical value of the bifurcation parameter κ 1 (or κ 2 ). Therefore, its trivial solution is a centre and its general steady state solution is where the overbar indicates complex conjugate and ω is given by Eq. (25) for a specific τ value and for κ 1 = κ 1,cr (or κ 2 = κ 2,cr ). c = [1 c 2 c 3 c 4 ] is obtained by the linear homogeneous system of equation where I is the identity matrix, which gives and for the collocated and non-collocated cases, respectively.
Because the homogeneous part of Eq. (40) has nontrivial solutions, the non-homogeneous equation has a solution only if a solvability condition is satisfied. To determine this solvability condition, we assume a particular solution of the form and obtain The solvability condition demands that the right-hand side of Eq. (42) is orthogonal to every solution of the adjoint homogeneous problem. The adjoint problem is governed by where H indicates the conjugate transpose. To make b unique, we impose the condition for the collocated and non-collocated cases, respectively. We take the inner product of b H with the righthand side of Eq. (42), i.e.
from which we get the expression

Normal form
We introduce the polar form where a and β are real functions of T 1 . Substituting Eq. (49) into (48) and separating real and imaginary parts, we get where Λ 1R , Λ 2R , Λ 1I and Λ 2I are the real and imaginary parts of Λ 1 and Λ 2 , respectively. Equation (50) provides information about the oscillation amplitude of the system in the vicinity of the bifurcation. In order to identify periodic solutions, we impose a = 0 and obtain the solutions where for the non-collocated case.

Bifurcation diagrams
In this section, we analyse the results obtained through the bifurcation analysis performed with the multiple scales method in the previous section. These analytical results are compared with direct numerical integrations and with bifurcation diagrams obtained through the NDDE-Cont continuation software for neutral differential equations [6], which is based on the continuation toolbox DDE-BIFTOOL for MATLAB [16]. Initially, results obtained with a classical linear controller (κ nl,1 = 0 or κ nl,2 = 0) are presented. Later, possible improvements in the bifurcation behaviour related to introducing a nonlinear term in the controller are studied. For the sake of brevity, the values of χ , r and γ are fixed at 0.05, 1 and 1, respectively, throughout the whole analysis performed in this section. Nevertheless, extending the obtained results to other parameter values is a relatively simple task left for the interested reader.

Linear control: collocated system
Equation (53) enables us to characterise the types of bifurcations occurring at the loss of stability. At this stage, we assume that the controller is purely linear (κ nl,1 = 0). Based on the calculated value of Δ, Fig. 5a illustrates that for μ = 1, which corresponds to a hardening nonlinearity of the stiffness, bifurcations are always subcritical in the analysed parameter range. Indeed, in Fig. 5b it can be recognised that, for the stability boundary corresponding to κ 1,cr > 1 (κ 1,cr < 1), Δ is positive (negative), which indicates that all corresponding Andronov-Hopf bifurcations are subcritical. The practical consequence of a subcritical bifurcation is that a branch of unstable periodic solutions exists within the stable region of the trivial solution. Generally, this limits the robustness of the stable trivial solution, causing the system to undergo large oscillations if subject to sufficiently large perturbations, even within the stable region of the trivial solution.
It can be noticed that the sign of Δ directly depends on the sign of μ. Therefore, considering a softening stiffness nonlinearity of the system (μ = −1), we expect to have supercritical bifurcations instead of subcritical. This is confirmed in Fig. 5c, d, which shows exactly the same features of Fig. 5a, b, but with reversed sign of Δ, which indicates that all bifurcations are supercritical. This feature could be easily predicted by looking at Eq. (53) for κ nl,1 = 0.
In Fig. 6, bifurcation diagrams obtained analytically (black lines) overlap with bifurcation diagrams obtained through the NDDE-Cont continuation software (blue and orange lines). The two procedures provide practically identical results at low amplitudes, which confirms the validity of the analytical procedure and that a softening nonlinearity produces subcritical bifurcations.
Since the analytical procedure neglects terms of order higher than ε 3/2 , it provides a bifurcation curve that is restricted to a parabola, making it impossible to detect any other trend of the branches of periodic solutions arising at the loss of stability. However, more complicated shapes of the branches can be easily tracked by extending the numerical continuation to higher values. Figure 6c illustrates that also in the case of softening nonlinearity (blue line), the branch of stable periodic solutions encounters a fold point at κ 1 ≈ 0.5, losing stability and turning back towards the stable region. This behaviour generates unstable periodic solutions within the stable region of the trivial solution, even in the case of supercritical Andronov-Hopf bifurcations. Nevertheless, this unstable branch has a very high amplitude, especially if compared with the red curve in Fig. 6c referring to the hardening case. This scenario suggests that the trivial solution, although not globally stable, is probably much more robust in the case of softening rather than hardening nonlinearity.
In this respect, Fig. 7 illustrates that, for fixed parameter values belonging to the stable region of the trivial solution, the system can either converge towards the trivial solution or not, depending on the initial conditions. In the case of hardening nonlinearity (μ = 1), Fig. 7a shows that selecting initial conditions y(t) = [0.735, 0, 0, 0] or y(t) = [0.740, 0, 0, 0] with t ∈ (−τ, 0], the system either converges towards the trivial solution or towards a much larger attractor (seemingly quasiperiodic). A somehow similar scenario occurs in the case of softening nonlinearity (μ = −1), as depicted in Fig. 7b for the same parameter values (except μ, of course). However, in this case, the system still converges towards the trivial solution for a significantly larger initial condition of x 1 . Additionally, the non-convergent time series is not attracted to any bounded solution, but it diverges unboundedly. This behaviour is probably related to a statistical loss of stability due to the softening stiffness, which becomes negative for excessively large |x 1 − x 2 | values and not to the unstable periodic solution detected by the numerical continuation. Clearly, for such large |x 1 − x 2 | values, the mathematical model adopted fails to describe the mechanical system under study. Indeed, the unstable periodic solution encompasses points such that (x 1 − x 2 ) 2 > γ 2 , for which the stiffness is negative, being therefore out of the region of validity of the mathematical model. Although rigorously proving that for hardening nonlinearity the trivial solution is less robust than for softening nonlinearity is very demanding because the system is infinite-dimensional, the performed analysis hints that a similar conclusion is plausible.
An interesting feature of the bifurcation diagram in Fig. 6c, for μ = 1, is that the branch of unstable periodic solutions reaches values of κ 1 close to zero (the fold occurs at κ 1 = 0.023), meaning that large amplitude vibrations might unexpectedly occur also for κ 1 far from the stability boundary.
From the stability analysis of Sect. 3, we deduce that for κ 1 < −1 an infinite number of eigenvalues has positive real parts. This occurrence makes it very hard for the numerical solver to track a periodic solution and perform a continuation in the vicinity of this boundary. For this reason, bifurcation diagrams for negative κ 1 values are depicted in Fig. 6b, d for τ = 0.8, where κ 1,cr is equal to −0.87 and relatively far from −1 (see, for instance, Fig. 5a). Apart from confirming the analytical results, i.e. that the bifurcations are subcritical for hardening stiffness and supercritical for softening, the bifurcation diagram of Fig. 6d illustrates that the branch of unstable periodic solution exists for a wide range of κ 1 negative values (up to κ 1 < −0.07). This observation implies that also for κ 1 < 0, unexpected large amplitude oscillations might occur for almost any κ 1 value in the case of hardening stiffness. Figure 6e illustrates the bifurcation diagram for the same bifurcation point of Fig. 6d, but adopting the time delay τ as bifurcation parameter. For the hardening case , the figure shows that the branch of unstable periodic solutions is very steep, rapidly reaching high values. This suggests that in this case, for low values of time delay τ , although the basin of attraction of the trivial solution is probably bounded, the trivial solution might still be globally stable limiting initial conditions to states of the mechanical system compatible with a physical system. The steep trend of the branch of periodic solutions for negative κ 1 could be already predicted by analysing the absolute value of Δ in Fig. 5b, d, which is tiny for small τ values. In fact, Eq. (55) clearly shows that the smaller the absolute value of Δ, the steeper the branch of periodic solutions arising at the bifurcation.

Linear control: non-collocated system
We now consider the non-collocated system and repeat the analysis performed in the previous section for the collocated system. Analysing the value of Δ from Eq. (54), the Andronov-Hopf bifurcations occurring at the loss of stability are characterised. Figure 8a, c illustrates that, similarly to the collocated case, also for the non-collocated system, the bifurcations are subcritical in the case of hardening stiffness and supercritical in the case of softening stiffness. Fixing the time delay τ at 0.15, bifurcation diagrams are depicted in Fig. 9. Also in this case, the good matching between analytical and numerical results exhibited in Fig. 9a, b validates the analytical procedure. Figure 9c, d provides a wider view of the bifurcation diagrams, which enables us to draw some conclusions about the effects of the bifurcations on the global system dynamics. In Fig. 9c, it can be noticed that both branches, for μ = 1 and μ = −1, increase their amplitude very slowly, if compared with the bifurcation diagrams in Fig. 6b. This observation could be deduced by the analytical procedure by looking at the absolute value of Δ in Fig. 8b, d. For κ 2,cr < 0 and decreasing τ , the absolute value of Δ rapidly increases, such that lim τ →0 Δ = ∞. This is exactly the opposite trend of Δ in the collocated case for κ 1,cr < 0, where Δ tends to 0 with a horizontal tangent for τ → 0 (Fig. 5b). Such a large Δ value, in the case of a supercritical bifurcation, makes arising periodic solutions have a relatively small amplitude (see blue curve in Fig. 9c), which is a x 1max x 1max x 1 x 1 t t desirable feature from an engineering point of view in many applications. On the contrary, if the bifurcation is subcritical, a low branch of unstable periodic solutions is generated, which might significantly reduce the robustness of the stable trivial solution. In Fig. 10a, time series obtained from direct numerical simulations illustrate the phenomenon. For μ = 1, κ 2 = −3 and τ = 0.15, therefore inside the stable region, a small perturbation of the value of x 1 (x 1 (0) = 0.27 in the figure) is sufficient to make the system diverge from the trivial solution causing unbounded vibrations. For the same parameter values but with a softening stiffness (μ = −1), the system converges to the trivial solution even for much larger initial displacement x 1 , as shown in Fig. 10b.
Considering bifurcations occurring at the positive stability boundary (κ 1,cr > 0), bifurcation diagrams for τ = 0.15 are illustrated in Fig. 9d. For the softening case (μ = −1), the branch has a fold, similarly to the bifurcation diagram encountered in the collocated case (Fig. 6a). However, in this case, the fold is very far from the stability boundary and does not seem to endanger the robustness of the trivial solution inside its stable region. Conversely, for a hardening stiffness (μ = 1), the branch of unstable periodic solutions is rather steep and does not enter very deeply in the stable region (excluding the part of the branch having huge amplitude), meaning that the region of limited robustness of the trivial solution might be not very large.
The time series in Fig. 10c proves that, in the vicinity of the stability boundary (κ 2 = 0.3, τ = 0.15) and for softening nonlinearity, the system is bistable. Namely, depending on the initial conditions, the system's dynamics either converges to the trivial solution (blue curve) or to a larger periodic solution (red curve), although the existence of other attractors cannot be excluded. We also remark that the amplitude of the large attractor (x 1max = 1.6) and of the unstable periodic solution on which the system initially lies (x 1max ≈ 0.52) in Fig. 10c agree with the results provided by the numerical continuation illustrated in Fig. 9d.

Bifurcation manipulation via nonlinear control
All bifurcation diagrams illustrated in the previous sections refer to the case of a linear controller, i.e. with κ nl,1 = 0 (and κ nl,2 = 0). However, Eq. (53) suggests that a proper choice of κ nl,1 might change the sign of Δ and switch the bifurcation characteristic from subcritical to supercritical or the other way around.
Δ is linear with respect to κ nl,1 , therefore values of κ nl,1 such that Δ is either positive or negative generally exist. Furthermore, the value of κ nl,1 causing the transition from subcritical to supercritical behaviour (i.e. such that Δ = 0) can be easily found analytically. Value of κ nl,1 and κ nl,2 at the bifurcation transition from subcritical to supercritical are illustrated in Fig. 11a, b for the collocated and non-collocated cases, respectively. Since bifurcations are already supercritical for softening stiffness, only the case of hardening stiffness (μ = 1) is considered, where forcing subcritical bifurcations to be supercritical might be beneficial.
The black curve in Fig. 11a refers to the upper stability boundary, while the red curve to the lower one. Interestingly, we notice that for the upper stability boundary a negative value of κ nl,1 is required, while for the lower one, a positive value. From an engineering point of view, this is not a desired feature, since, for instance, a positive value of κ nl,1 might be beneficial from the point of view of the global stability for κ 1 > 0, but it will be detrimental for κ 1 < 0. Another aspect highlighted by Fig. 11a is that just a very small value of κ nl,1 is required to change the character of the bifurcation of the lower stability boundary, while for negative values of κ 1 larger absolute values of κ nl,1 are required. This is directly related to the Δ value depicted in Fig. 5b.
Regarding the collocated case, Fig. 12 illustrates the effect of the additional nonlinearity in the control law on the bifurcation diagram. Figure 12a refers to the upper stability boundary for τ = 0.1. According to Fig. 11a, the transition from subcritical to supercritical behaviour occurs at κ nl,1 = −0.0093. Figure 12a confirms this prediction, showing that for κ nl,1 = 5 × 10 −5 , 0, −5 × 10 −5 , −5 × 10 −4 , −0.005 the bifurcation is always subcritical, while it becomes supercritical for κ nl,1 = −0.01. The different bifurcation diagrams in the figure illustrate that negative values of κ nl,1 , in absolute value even smaller than 0.0093, generally reduce the extent of the bistable region, making the stable trivial solution more robust. Conversely, if κ nl,1 is positive, the bifurcation branch moves to the left. This effect most probably reduces the robustness of the stable region. Figure 12b depicts bifurcation diagrams occurring at the lower stability boundary, for κ 1 < 0 and τ = 0.8. According to Fig. 11a, the transition of the bifurcation from subcritical to supercritical should occur at κ nl,1 = 0.00033. This prediction is also confirmed numerically in the figure, where the bifurcation diagram for κ nl,1 = 0.0004 presents a supercritical characteristic. In this case, although the analytical procedure can correctly predict the type of bifurcation, branches of unstable periodic solutions exist within the stable region also in the case of supercritical bifurcations (as for κ nl,1 = 0.0004). Nevertheless, further increasing κ nl,1 (κ nl,1 = 0.001) these branches seem so disappear.
From Fig. 12a, b, it can be noticed that bifurcation branches can move from negative to positive ranges of κ 1 and vice-versa. This feature, which was not observed τ τ in the case of a linear control force, implies that the tuning of κ nl,1 should be done with particular care. In fact, if the linear control parameters are such that the system is, for instance, close to the upper stability boundary and we aim at enforcing supercritical bifurcation to grant robustness to the system, we need to tune κ nl,1 to a negative value. However, this will cause the bifurcation branch arising at the lower stability boundary to bend to the right, reaching parameter spaces where κ 1 is positive, potentially mining the system's robustness in the point of interest. This feature is better illustrated below for the non-collocated case. Figure 11b, referring to the non-collocated system, depicts the values of κ nl,2 such that the bifurcations occurring at the loss of stability undergo a transition from subcritical to supercritical. Comparing Fig. 11a, b, we notice that, for both the collocated and noncollocated cases, a positive nonlinear coefficient is required to force bifurcations at the lower stability to be supercritical (except for τ < 0.045) and a negative one for the bifurcations at the upper stability boundary. However, in the non-collocated case, much larger controller nonlinearity values are required, particularly for small τ values. This peculiarity is related to Δ, which reaches very large values for κ nl,2 = 0, as visible in Fig. 8b.
In Fig. 13, bifurcation diagrams with various values of κ nl,2 are represented for the non-collocated case. In the figure, τ is fixed at 0.15. According to the analytical prediction, the transition of the bifurcations from subcritical to supercritical should occur at κ nl,2 = 9.3 at the lower boundary and κ nl,2 = −0.69 at the upper boundary. This prediction is confirmed by the bifurcation diagrams in Fig. 13.
Apart from correctly estimating the bifurcation transition, Fig. 13a also shows that, in general, positive x 1max x 1max values of κ nl,2 move the branch of periodic solutions to the left, probably increasing the robustness of the trivial solution, even if κ nl,2 is smaller than the critical value 9.3. On the contrary, in the upper stable region (κ 2 > 0), because of the shape of the branch of periodic solutions, characterised by multiple folds, it is arguable if the additional nonlinearity produces any beneficial effect in terms of robustness. In Fig. 13b, red curves depict bifurcation branches arising at the lower stability boundary but reaching a region of positive κ 2 . Interestingly, these curves lie below the bifurcation branches arising at the upper stability boundary for the same κ nl,2 values, mining the robustness of the trivial solution. This result highlights that introducing an additional nonlinearity in the control force must be done with particular care and considering the system's global behaviour. This phe-nomenon is particularly relevant in the non-collocated case, where large nonlinearities are required to control bifurcations. Indeed, it is well known that systems encompassing large nonlinearities might present dynamical phenomena hardly predictable. We remark that bifurcation branches crossing the border of κ 1 = 0 (or κ 2 = 0) are present also in the other cases studied in this section but are not represented for clarity of illustration of the figures.
In order to directly illustrate the potentiality and limitations of utilising a nonlinear control force for enhancing robustness, Fig. 14 depicts time series obtained from direct numerical simulations. Figure 14a refers to the collocated system and parameter values close to the upper boundary (same parameter values as in Fig. 7a). If κ nl,1 = 0 and the initial displacement x 1 is too large, the system converges towards a large quasiperiodic attrac- tor (red curve). However, including a small negative nonlinearity in the control force (κ nl,1 = −0.005), for even larger initial conditions, the system still converges to the trivial solution. This result clearly suggests that the trivial solution is more robust thanks to the additional nonlinearity.
Referring to the non-collocated case, Fig. 14b shows time series for parameter values close to the upper stability boundary (same parameter values as in Fig. 10c). In this case, for certain initial conditions, the system either converges towards a large periodic attractor or towards the trivial solution, depending on the value of κ nl,2 (red and black curves in the figure, respectively, the large periodic attractor is visible in Fig. 10c). However, if the initial x 1 displacement is slightly increased, the system diverges unboundedly (blue curve) in the case of a nonlinear controller. This example directly shows that, although a nonlinear control force can increase robustness, it might also generate undesired and unexpected phenomena.

Conclusions
This study investigated the stability and bifurcation behaviour of a two-DoF system subject to delayed acceleration control. Two configurations of the same systems were considered, namely a collocated and a non-collocated one. The stability analysis illustrated the importance of considering the time delay of the control force, especially for the collocated case, where even an infinitesimal time delay can cause instabilities. Additionally, time delay seems particularly critical for both the collocated and non-collocated cases if damping is small.
The bifurcation analysis showed an interesting and unexpected feature about the global dynamics of the system. For both configurations studied, if the system encompasses a hardening stiffness nonlinearity, bifurcations at the loss of stability are subcritical, while they are supercritical in the case of softening nonlinearity. This behaviour has important consequences on the robustness of the stability of the trivial solution. In fact, particularly in the case of hardening nonlinearity, the system can experience large oscillations even within the stable region of the trivial solution if subject to a sufficiently large perturbation.
By adding a nonlinear term in the control force, a strategy to increase the robustness of the trivial solution was developed. Our analysis illustrates that, for a collocated controller, even a small nonlinearity can change the character of the bifurcations, improving the robustness of the trivial solution, even in the case of hardening nonlinearity. However, this approach seems to have significant limitations in the non-collocated case, where strong additional nonlinearities would be required to control bifurcations; these strong nonlinearities typically generate other undesired dynamical phenomena, hardly predictable.
The bifurcation analysis was performed for a given set of parameter values. The effect of the various parameters on the system's nonlinear behaviour might illus-trate other dynamical mechanisms not disclosed in this study. Other relevant aspects omitted in this work are the saturation of the control force and external excitations. These two elements might have an essential effect on the global dynamics of the system and will be analysed in future studies.
Funding Open access funding provided by Budapest University of Technology and Economics. This research work was done under the framework of the project MIRAGED: Posicionamiento estratégico en modelos virtuales y gemelos digitales para una industria 4.0 (CER-20191001). This research was supported by the InterQ European project (H2020-958357). GH acknowledges the financial support of the Hungarian National Science Foundation under grant number OTKA 134496. The research reported in this paper was supported by the National Research, Development and Innovation Fund (TKP2020 NC, Grant No. BME-NCS and TKP2021, Project no. BME-NVA-02) under the auspices of the Ministry for Innovation and Technology.

Data availability
The data sets generated during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors declare that they have no conflict of interest.
Code availability The code utilised for generating the presented results are available from the corresponding author on reasonable request.
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/.