Bifurcations in basic models of delayed force control

Basic single-degree-of-freedom mechanical models of force control are presented to achieve desired contact forces between actuators and objects. Nonlinear governing equations are constructed for both collocated and non-collocated force sensor configurations. The models take into account large time delays in the feedback loop, which may occur, for example, in case of human or remote force control. The corresponding stability charts are compared for the collocated and non-collocated cases. The bifurcations at the stability boundaries are analyzed in the presence of the relevant nonlinearity that is originated in the saturation of the actuation. The stability properties as well as the nonlinear vibrations for the two sensor locations are compared also from the viewpoint of the achievable maximal proportional gains. The bifurcation calculations are done with the method of multiple scales and normal form calculations. The results on global dynamic properties are supported with numerical simulations, and the two force control strategies are discussed from application viewpoint.


Introduction
The tasks of human motion control can be classified in many different ways. One of the simplest categorizations can be represented with the help of the basic single-degree-of-freedom (DoF) mechanical model where m and k refer to the inertia and the stiffness of the controlled mechanical system, respectively, which may be the model of certain parts of the human body, while Q refers to the control force of certain muscles, which is actuated as a function of the sensed position q, velocityq and/or accelerationq with dot referring to time derivatives. The position control task corresponds to the case when the mass is not in contact with the environment, that is, the stiffness is zero: k = 0. Force control is represented by the case of positive stiffness k > 0 that models the elastic contact between the mass m and the environment. Finally, stabilization of otherwise unstable positions refers to the case of k < 0, since these systems would be unstable without any control force Q due to the negative stiffness in the structure. Human balancing during standing still belongs to this class of problems.
Clearly, many further combinations of these basic control tasks can appear, like hybrid position/force control, or following prescribed paths while balancing, or tasks with many DoF involved. These basic models are established in classical robotics papers [1,2] and books like [3]. Biomechanics often uses these models for the mathematical description of human motion, although there are many obvious differences. One of these is related to the time delay that appears in the human motion control. The delays are originated in the processing time of the sensory signals in the neural system and also in the reaction time of the muscles. These time delays of 0.1-1 s are orders of magnitudes larger than the ones that appear in robotics related to the sampling times of the digital controllers. The digital effects related to sampling and zero-order hold in time were considered in several studies on force control [4][5][6][7], but it is less investigated in case of analogous time delays that are essential elements of models in human motion control [8], human force control [9,10] or human balancing [11][12][13][14]. Remote force control tasks of telemanipulation [15][16][17] and physical contacts between robots and humans via force control are actual research topics of haptic systems [18][19][20][21], specifically in rehabilitation robotics [22][23][24].
This paper investigates the basic models of force control in the presence of human reaction time delays from nonlinear vibrations viewpoint similarly to the Hopf bifurcation analyses in delayed oscillators [25][26][27], predator-prey population models [28][29][30] or delayed networks in transportation, biological and neural systems [31][32][33]. These calculations are usually based on the center manifold reduction and normal form theory summarized in [34], which are then combined with the infinite dimensional representation of delay differential equations [35]. The corresponding algorithm of the Hopf bifurcation analysis in delay systems was published in [36]. While this algorithm has a rigorously proven theoretical background, the application of the method of multiple scales (MMS) can make shortcuts in the lengthy algebraic manipulation [37]. A relevant single-authored paper [38] of Nayfeh discussed the comparison of the two approaches in the presence of delay parameters, and since then, MMS has become equally popular in the nonlinear analysis of delay systems.
The structure of the paper is as follows. The simplest possible force control strategy is considered with proportional gains while the saturation of the control force is taken into account as the relevant nonlinearity. The corresponding mechanical models are presented in Sect. 2, where the equations of motion are derived for different locations of the force sensors: first, called non-collocated case, when the contact force is sensed between the rigid environment and the mass, and second, called collocated case, when it is sensed between the mass and the actuator. In the non-collocated case, in Sect. 3, the stability analysis of the steady-state desired contact force is presented together with the Hopf bifurcation calculation done by the method of multiple scales extended for nonlinear delay differential equations. The same calculations are presented for the collocated configuration in Sect. 4. The analytical results and conclusions are summarized in Sect. 5.

Mechanical models of basic force control concepts
Consider the basic task of force control when a block of mass m is in contact with the rigid environment through a spring of stiffness k 1 along a horizontal axis as shown in Fig. 1. The actual contact force is measured by means of a serially connected spring of large stiffness k 2 ( k 1 ). In Fig. 1a, the so-called non-collocated configuration is presented where the force sensor is at the contact point to the rigid environment and its signal is fed back to the control force Q of the actuator, while Fig. 1b is the sketch of the collocated configuration with contact force measured at the actuator force.

Modeling the non-collocated configuration
In case of Fig. 1a, the Newtonian equations assume the form where the selected generalized coordinates q 1 and q 2 are the absolute positions of the block and the end point of the spring that senses the force, respectively. The actuator force Q depends on the force error Mechanical models with non-collocated (a) and collocated (b) force sensor configurations. The stiffness k 2 of the force sensor is an order of magnitude larger than the stiffness k 1 of the actuator. P refers to proportional gain; τ stands for time delay which is the difference in the actual measured force and the constant desired force F d . Consider the simplest possible linear control strategy using a single (dimensionless) proportional gain P for the force error F e while adding the actual measured force (see [3]): The reaction time τ of the human actuation can be taken into account with the delayed values F m (t − τ ) of the measured force in (4), and the additional parameter F s can describe the level of the actuator force saturation approximated by means of the tangent hyperbolic function tanh. This way, the actual control force Q at the time instant t is given by The equations of motion (2) with the nonlinear delayed control force (6) have an equilibrium at the trivial solution The perturbation x around the equilibrium of the block at q 10 with q 1 (t) = q 10 + x(t) is introduced. If the approximation k 2 k 1 is used for the angular natural frequency ω n of the uncontrolled system, and if the dimensionless perturbed block positionx is scaled to the static saturation level F s /k 1 and the dimensionless timet is scaled to the delay τ , then we have where, by abuse of notation, we will drop the tilde from the dimensionless variables. Substituting all these into the Newtonian equation (2) and eliminating the state variable q 2 , we obtain the nonlinear delay differential equation (DDE) of retarded type in the dimensionless form of with trivial solution at zero. Note that the dimensionless parameter ω n τ = 2π(τ/T ) can also be expressed with the ratio of the reaction time τ and the time period T of the free oscillation of the uncontrolled system.

Modeling the collocated configuration
When the collocated configuration is considered as shown in Fig. 1b, the Newtonian equations assume the form where the generalized coordinates q 1 and q 2 now stand for the absolute position of the block of mass m and the position of the end of the spring of stiffness k 2 relative to the position of the block. This choice of coordinates means that the measured force F m can be expressed in the same way as in (4) and the saturated and delayed actuator force can be expressed as given in (6). The equations of motion (10) with the control force (6) have a trivial solution at If the perturbation x is introduced again around the equilibrium q 10 of the block in the same way as in the case of the non-collocated configuration, and the same dimensionless quantities are used as in (8), then the elimination of the coordinate q 2 from the equations of motion (10) leads to the nonlinear DDE of neutral typë This equation of motion is similar to Eq. (9) of the non-collocated case, but there is an essential difference due to the appearance of the delayed acceleration terms x(t − 1). While (9) has the standard form of a retarded DDE, the neutral DDE (12) has an implicit form that generally cannot be expressed in the standard explicit form as required for the application of fundamental theorems for neutral DDEs (see, e.g., in [35,39]). In this case, however, it is easy to recognize that the implicit neutral DDE (12) can be transformed into the explicit scalar nonlinear difference equation by means of the new variable The physical meaning of this transformation is clear: The dynamics of the block and the spring of stiffness k 1 has no effect on the dynamics of the actuator force that is determined by the delayed force k 2 q 2 (t − τ ) sensed at the actuator. In other words, the second algebraic equation of (10) with the delayed and saturated control force expression in (6) can be transformed into the same nonlinear difference Eq. (13) by selecting y = k 2 q 2 − F d .

Analysis of non-collocated configuration
The third-order Taylor series approximation of the nonlinear retarded DDE (9) assumes the form First, the linear part is used to determine the stable domains in the parameter plane ω n τ and P, and then, the Hopf bifurcation calculations are presented for the above third-order approximation by means of the methods of multiple scales.

Stability
The local stability is analyzed by means of the characteristic function of the linear part of (15). At the limit of stability, the decomposition of the characteristic function (16) with substitution of the critical pure imaginary characteristic root λ c = iω c yields: The equations Re(D(iω c )) = 0 and Im(D(iω c )) = 0 give critical values for the gain P and the dimensionless critical angular frequency ω c : with k = 1, 2, . . . . This leads to the stability chart in the parameter plane of (P, ω n τ ) for the case of noncollocated force control, as shown in Fig. 2. For certain given values of the gain P, the critical delay parameters τ c can be obtained from where k = 1, 2, . . ., and the critical angular oscillation frequency is given by ω c, j = jπ corresponding to the stability limit at τ c, j . Some numerical values for the critical delays τ c, j and vibration frequencies ω c, j are summarized in Table 1 for the specific gains P = 0.6, 0.9, 1.2 and 1.5 in the subsequent subsection where the Hopf bifurcations are analyzed at the critical delays (20) and (21). The stability boundaries in Fig. 2 are numbered according to the sequential numbers j of the critical delays, but for the sake of simplicity,  As it is shown in Fig. 2, the maximal gain value can be calculated at the intersection of the first and second stability boundary curves. The calculation is based on the second and third equations of (19) for k = 1, and it gives P max = 8/5 at the critical delay τ c,1 = τ c,2 = √ 5/2 π/ω n .

Hopf bifurcation
In order to analyze the effect of varying the delay parameter, the dimensionless bifurcation parameter μ is chosen as in the subsequent calculations, where the critical dimensionless parameter values ω n τ c corresponding to the critical delay τ c can be determined from (19) for any fixed value of the gain P. To study small amplitude oscillations, let where ε is a small non-dimensional bookkeeping parameter, i.e., 0 < ε 1, andμ = O(1) is the detuning parameter (see [37,38] or Appendix 2 in [40]). Then Eq. (15) can be transformed in the form where, by abuse of notation, the over-lines are dropped immediately from x and μ. The multiple timescales are defined as T k = ε k t, k = 0, 1, 2, . . .. To study the Hopf bifurcation, a two-scale expansion of the solution is assumed as By means of the differential operators defined as [38] d the delayed state can be expressed as Substituting Eqs. (25), (26) and (27) into Eq. (24), and equating the same powers of ε, a set of linear partial differential equations are obtained in the form where the notations x k := x k (T 0 , T 1 ), x kτ := x k (T 0 − 1, T 1 ) are introduced. At the stability boundary, one pair of pure imaginary characteristic roots ±iω c exists, while all the other infinitely many eigenvalues have negative real parts. All the solution terms related to these negative real eigenvalues decay with time. Thus, to study the long-time behavior of the system, the solution of the first equation of (28) is assumed as with the complex amplitude R and with R denoting its complex conjugate, i 2 = − 1. Substituting Eq. (29) into the second equation of (28), the secular term can be found, which must be zero: where over-circle stands for the derivative w.r.t. the slow time: This can be arranged into the complex normal form The sign of the Poincaré-Lyapunov constant Re δ determines the sense of Hopf bifurcation: It is supercritical or subcritical if Re δ < 0 or Re δ > 0, respectively, i.e., the bifurcated periodic motions are stable or unstable, respectively. Here for the non-collocated case, we have which indicates that the Hopf bifurcation is supercritical (or subcritical) for P > 1 (or 0 < P < 1). These results are also indicated in case of the numerical values in Table 1. Note that there are degenerate Hopf bifurcations along the stability boundary P = 1 at the Bautin points ω n τ = jπ , j = 1, 2, . . ., where not only the sense of the Hopf bifurcations, but also the stability of the equilibrium change to opposite in the same time. In this respect, these Bautin points are special ones.
The bifurcating periodic motions are approximated by the dimensionless form that can be obtained from (29) when the phase is eliminated by the appropriate selection of the initial time instant. The dimensionless vibration amplitude r = 2|R| and the dimensionless angular frequency ω c satisfy (31) with • R ≡ 0 and R = 0, that is, Its real and imaginary parts lead to the equations sin ω c = 0.
The second Eq. (38) is satisfied at the stability boundaries given in (19) by ω c = jπ , ( j = 1, 2, . . .). The substitution of cos ω c = (−1) j into (37) gives |R|, the double of which gives the dimensionless vibration amplitude in the approximation of the periodic solutions (35). This means that the branches bifurcating from the first, third, fifth, . . . stability boundary curves (for j = 2k−1) have amplitudes while for the bifurcating branches from the second, fourth, sixth, . . . boundary curves (for j = 2k) are defined by and they are all unstable for 0 < P < 1 and stable for 1 < P < 2. The corresponding bifurcated branches are presented in Fig. 3 by means of the above algebraic formulas, and these are also checked by numerical solutions of the original nonlinear Eq. (9) obtained with a collocation method [41] combined with path following.
Transforming these results back to the time and space with dimensions according to (8), approximation (35) gives the periodic contact force oscillations bifurcating from the jth stability boundary in the form: with j = 1, 2, . . .; these are all unstable for 0 < P < 1 and stable for 1 < P < 2. In this formula, the jth critical delays τ c, j are given in (20), (21) and the delay parameter τ is 'close enough' to one of these. Recall that in formula (42), F d is the desired contact force and F s is the force saturation level of the actuator. The comparison of the numerical and analytical results in Fig. 3 shows that the algebraic approximations obtained by MMS are acceptable in engineering applications for delay parameters τ within 15% of their corresponding critical values. The Hopf bifurcations can also be studied with the same methodology along the stability boundary P = 1 given in (19), too. In this case, the modified definition of the bifurcation parameter μ is necessary along the gain parameter P in the form Since the linear system does not include any delay term in the critical parameter case μ = 0, the calculations with MMS are simpler than the ones above. For the sign of the Poincaré-Lyapunov constant, one gets Sgn Re δ = Sgn (− sin(ω n τ )).

Analysis of collocated configuration
When the collocated force control configuration is analyzed, the relevant difference Eq. (13) is considered with the third-order polynomial approximation Since this is a simple scalar nonlinear difference equation, the stability and bifurcations of the zero solution can easily be obtained based on the normal form of scalar discrete maps [34].

Stability
The stability of the linear part of (45) is equivalent to the convergence of the geometric series of quotient (1− P), which is the characteristic multiplier of the corresponding scalar linear map. The stability condition assumes the form which means that the stability of the collocated configuration does not depend on the delay, and the maximal proportional gain that can be achieved without losing stability is P max,c = P c = 2 as opposed to the non-collocated configuration where the maximal gain P max,nc = 1.6 can be achieved at τ = √ 5/2 π/ω n only.

Period-doubling bifurcation
At the critical gain P c , the characteristic multiplier is − 1, that is, the nonlinear mapping (45) has a perioddoubling (or flip) bifurcation here. Since only the gain parameter is responsible for the loss of stability in this case, the bifurcation parameter μ is now defined by and the normal form of the corresponding mapping is given as Since δ = 8/3 > 0, the period-doubling bifurcation is supercritical, that is, stable period-2 oscillations appear in (45) and similarly in (13) for μ > 0, that is, for P > P c . The amplitude r of this oscillation is expressed by After transforming this result back to the real contact force in the real time, it means that for proportional gains P somewhat larger than the critical value P c , stable self-excited contact force oscillation occurs around the desired value with frequency 1/(2τ ) that corresponds to the basic critical dimensionless angular frequency ω c = π related to the critical characteristic multiplier at exp(iω c ) = −1:

Conclusions
Basic single-degree-of-freedom mechanical models of force control with collocated and non-collocated force sensor arrangements were constructed and the governing delay differential equations (DDE) were derived when the relevant nonlinearity is the actuator saturation and time delay appears in the feedback loop. These models are similar to those of position control with collocated and non-collocated position sensors [42], but the equations of motion are more complex due to the unilateral elastic connection of the controlled object and the environment. The mathematical model of the non-collocated force control is a retarded DDE that belongs to the class of delayed oscillators. The stability properties show alternating stable and unstable regions for both the delay and the gain parameter (see Fig. 2). The maximization of the gain P in the feedback loop within the stable parameter region has relevance when there is dry friction in the mechanical structure. The maximum gain P max = 1.6 can be achieved only at the critical delay τ c = √ 5/2 π/ω n where ω n is the angular natural frequency of the uncontrolled system. At the stability boundaries, subcritical Hopf bifurcations occur for proportional gains P < 1, while these are supercritical for P > 1. At loss of stability, the vibration frequency and the vibration amplitude of the contact force are well approximated with the simple algebraic formula (42).
The mathematical model of the collocated force control is an implicit neutral DDE that can be transformed into a simple nonlinear difference equation. Its dynamic properties are much simpler: The system is stable for all gains 0 < P < 2, and the Hopf bifurcation is always supercritical at the gain P max = 2. This maximal gain can be achieved for any time delay, which means that the system is robust against the delay parameter. In practice, however, the accuracy of this control strategy is worse due to the fact that the sensor is not exactly at the position of the contact between the object and the environment. The approximation of the emerging stable oscillations for P > 2 is presented again by the simple formula (50).
The closed-form analytical expressions for the stable parameter domains are useful for the preliminary parameter estimation during the design of force controlled structures. The formulas for the frequencies and amplitudes of emerging vibrations support the identification of experimentally detected unwanted vibrations and also the parameter estimations regarding the critical ratios of the time delays and the time periods of the natural oscillations of the uncontrolled system.
There are several open questions even in these simple models of force control, especially in the noncollocated case. The numerical results show that the bifurcated branches leaning against each other between subsequent stability boundaries remain separate as it is shown, for example, by the rightmost unstable branches in Fig. 3 for P = 0.9. These periodic solutions are probably 'separated' by the insets and outsets of quasiperiodic oscillations embedded in the infinite-dimensional phase space. The structure of these could be explored with the methods presented in [34,43,44] for the analysis of codimension-2 Hopf bifurcations that emerge at the cross sections of subsequent Hopf stability boundaries shown also in Fig. 2. Among these, the one at the peak of the stability chart in Fig. 2, where the gain P is maximal, requires even more careful handling since this codimension-2 Hopf bifurcation occurs with the specific 1:2 resonant frequencies.
The gains 0 < P < 1 seem to guarantee exponential stability for very short time delays and/or for very soft connections between the controlled object and the environment. However, the subcriticality of the Hopf bifurcation at the first bifurcation point for the delay shows that the force control can still exhibit increasing vibrations for large enough perturbation. There are no attractors identified outside these unstable branches, and no large-amplitude stable oscillations are identified for increased values of the delay either (see the bifurcation diagrams for P = 0.6 and 0.9 in Fig. 3). In these cases, another relevant nonlinearity has to be included and analyzed in the model: The increasing large-amplitude oscillations are limited by the loss of contact between the controlled object and the environment. Such large-scale nonlinearities originated in the unilateral constraints are relevant in cases of transitions between position and force control (see, e.g., [45,46]), but the study of these global nonlinearities in the noncollocated force control model for 0 < P < 1 is the task of future research.
Another possible extension of the research work is the study of the effect of viscous damping in the model. It is expected that the disjunct stable domains of the non-collocated case become simply connected, the stable region gets larger via a smooth variation of the stability boundaries with increasing damping, the special Bautin points disappear, while the double Hopf bifurcation points still survive. In the collocated case, the separated oscillator becomes physically damped, too. However, the algebraic calculations become so intricate that the effect of the parameters on the nonlinear system behavior cannot be summarized in simple closed formulas anymore.