Abstract

Precision and stability in position control of robots are critical parameters in many industrial applications where high accuracy is needed. It is well known that digital effect is destabilizing and can cause instabilities. In this paper, we analyze a single DoF model of a robotic arm and we present the stability limits in the parameter space of the control gains. Furthermore we introduce a nonlinearity relative to the saturation of the control force in the model, reduce the dynamics of the nonlinear map to its local center manifold, study the bifurcation along the stability border and identify conditions under which a supercritical or subcritical bifurcation occurs. The obtained results explain some of the typical instabilities occurring in industrial applications. We verify the obtained results through numerical simulations.

Introduction

Position and force control in robotic systems have been the subject of intensive investigations in recent years. Robotic systems are usually equipped with digital controllers, while their dynamic analysis is often treated using continuous-time (analog) approaches and models [1–3]; indeed, the sampling frequency can be so high, that from an engineering point of view it can be considered continuous. However, these approaches do not provide always conservative estimates for stability.

In the last decade several authors studied the stability of digital control systems, using mainly linear models [4–7], finding surprising results about the stability domains, which are strongly influenced by the sampling frequency.

Industrial robots often present instabilities that cause the arm to oscillate around the desired position. The stability is then reached through several tests aimed at ‘‘calibrating’’ the control system. Even though considerable work has been devoted to analyzing bifurcations of discrete-time systems [8–10], up to the authors’ knowledge a detailed analysis of the bifurcations happening at the border of the relevant stability domains is still incomplete and requires further investigations.

In this work, we analyze the stability of the digital position control of a single DoF system, using a linear model and considering an infinitely rigid connection element between the actuator and the manipulator of the arm. The control system is a proportional-differential (PD) one. Then we introduce in the model a nonlinearity relative to the saturation of the control force, reduce the dynamics of the nonlinear map to its local center manifold, study the bifurcation along the stability border and identify conditions under which a supercritical or subcritical bifurcation occurs.

In the second part of the work, focusing on the case of asymmetric control force, we perform systematic numerical simulations of the system aimed at validating the results obtained analytically and at obtaining some information on the global response of the system. We define the region in the parameter space where the bifurcation is supercritical or subcritical, highlighting the coexistence of attractors, and discuss the onset of chaotic motions. We verify that the motion is actually chaotic based on the Lyapunov exponent and the fractal dimension of the attractor. The paper ends with some conclusions and suggestions for future work.

Position Control

Continuous Model.

Figure 1 presents an ideal one DoF mechanical model of the position control, where m stands for the mass modeling the inertia of the robot and Q is the control force. The desired position xd is set to zero in order to simplify the equations without losing generality of the problem. The digital processor reads the input, the position x and the velocity x· of the mass, and elaborates a current signal that is sent to the DC motor.

The characteristic curve of the force of the DC motor is considered linearly proportional to the input current, so we can express the control force as
(1)
where P and D are the proportional and the differential gains, respectively. The mass is considered to be rigidly connected to the motor. The equation of motion of the system is then
(2)

Clearly the trivial solution x=0 satisfies Eq. (2). Neglecting the dry friction C, this position is asymptotically stable if and only if P>0 and D>0. As shown in a continuous time approach, the stability domain is then unlimited and there is no need for further analysis.

Discrete-Time Model.

The digital processor samples the force signal at the time instants tj=jτ,j=0,1,2,..., where τ stands for the sampling time, f=1/τ is the sampling frequency. The control force is piecewise constant in the sampling intervals (this effect is also referred to as zero-order hold) and it is calculated from the position and the velocity sampled at the beginning of the previous sampling interval:
(3)
Neglecting the dry friction again, the equation of motion is
(4)
We then introduce the dimensionless time T as
so the equation of motion becomes
(5)
where prime stands for the dimensionless time derivative. Calling p=Pτ2/m and d=Dτ/m we can rewrite Eq. (5) in the simpler form
(6)
where the only parameters controlling the dynamics of the system are p and d. We define aj as an acceleration
(7)

which is piecewise constant.

Stability Analysis

Discrete Mapping.

After integration of equation, Eq. (7), we can define the acceleration, the velocity and the position in the following way:
(8)
(9)
(10)
where T-j, calculated in j + 1, is equal to 1. So in matrix form we have
(11)

Equation (11) is a discrete map that describes the dynamics of the system; clearly, the convergence of the 3D geometrical series, Eq. (11), is equivalent to the asymptotic stability of the zero solution of Eq. (4).

Convergence.

As explained in basic textbooks, the trivial solution of the linear mapping, Eq. (11), is asymptotically stable if and only if all the eigenvalues of A, i.e., the eigenvalues |μ1,2,3|<1 are located within the unit circle of the complex plane. However, the study of the possible bifurcations and vibration frequencies along the stability limits also requires the use of the characteristic exponents s. Their relations are represented in Fig. 2 in accordance with the following standard calculation [6].

The substitution of the exponential trial solution x(T)=KesT(sC,KR3) into Eq. (11) yields the characteristic equation det(esτI-A)=0. This has an infinite number of roots sk,k=1,2,... called characteristic exponents, which are situated along a finite number of vertical lines in the complex plane as shown in Fig. 2. Clearly, there is only a finite number of characteristic multipliers defined by μ=esτ. These are actually the three eigenvalues of the transition matrix A. The criterion for the exponential stability of the force control can be written as
(12)

The same stability condition can be derived using the Laplace transformation L, which is illustrated in Fig. 2.

For the stability analysis, we follow the calculation of the eigenvalues μ from the characteristic Eq. (12)
(13)
where
In order to avoid the tedious algebraic work during the solution of the above third degree equation, the so-called Moebius transformation can be used for the new variable σ [6,7].
(14)
The transformation is also shown in Fig. 2; the open unit circle is mapped back to the left half of the complex plane. After the application of the transformation in Eq. (14), the equation again becomes a third degree polynomial
(15)
(16)

the kth column of the transformation matrix contains the binomial coefficients of (σ+1)4-k×(σ-1)k-1, k=1,2,3,4.

In order to verify when the roots of p3 are located on the left half of the complex plane, the Routh-Hurwitz criterion can be applied to the coefficients in Eq. (16).

Stability Charts.

The Routh-Hurwitz criterion gives the following conditions to be satisfied in order to have stability:

In Fig. 3, the stability chart is shown. The maximal allowed value of p, in order to have stability, is p=1/4 when d=5/8, i.e., Pτ2/m=1/4 and Dτ/m=5/8. From the stability chart we immediately notice that the stability domain is bounded due to the digital effect.

Dynamic Behavior

The most interesting stability limit in Fig. 3 is the one where H2=0, since the other (p=0) is a simple static loss of stability. It is easy to find that when the coefficient p is increased up to the stability limit, two complex conjugate eigenvalues of A go out from the unit circle of the complex plane with the value μ1,2=e±iα. According to Floquet theory in this circumstance there is a Neimark-Sacker bifurcation [11].

Vibration Frequencies.

In the case of the Neimark-Sacker bifurcation, the lowest angular frequencies of the oscillation at the loss of stability can be calculated in closed form. To do so we substitute the generic value iω to σ in Eq. (15), then we separate the real and the imaginary parts in order to find the value of ω on the curve H2=0. When σ1 has a pure imaginary value, we can easily find that
(17)
Since |μ1|=1 and μ1=|μ1|eiα=eiα=cosα+isinα, so
(18)
The vibration frequency is then
(19)

From Eq. (19) we plot the vibration frequencies along the stability limit in Fig. 3.

Nonlinear Analysis

Nonlinear System.

In order to study the bifurcation at the stability limit, we consider the nonlinearity arising from the saturation of the control force in the equation of motion. We approximate the force of the DC motor as an arctangent function of the current, while the current is a linear function of the displacement and the velocity of the mass m. Eq. (1) becomes
(20)

where Qmax>0 is the maximal output force of the motor, the parameter q0 indicates that the saturation of the control force is not symmetric, while the parameter u0=tanq0 is introduced in order to have zero force for u=0 (Fig 4).

We apply the transformation T=t/τ and we call Q0=2Qmax/(mπ), so we can rewrite Eq. (20) as
(21)
Transforming the arctangent in a Taylor series we have
(22)
and collecting terms with the same power order
(23)
In order to have the same stability domain as the one in Fig. 3 and the same linear part, we redefine p=τ2PQ0/(1+u02) and d=τDQ0/(1+u02) the coefficients relative to the linear part, while we call pij the coefficients relative to the nonlinear part, in the following way:
(24)
So we have
(25)
(26)
(27)
and the nonlinear map
(28)
where
(29)

Jordan Normal Form.

In order to study the bifurcation on the stability limit, we focus on the point with the maximum allowed value of p. So we fix d=5/8 and we study the bifurcation for the variation of p; we define H=A|d=58. The eigenvalues of H for ppcr=1/4 are
(30)
The corresponding eigenvectors are
(31)
Using the eigenvectors in Eq. (31) we define the transformation matrix
(32)
from which we can easily verify that
Applying the transformation
(33)
we can rewrite system, Eq. (28), in Jordan normal form
(34)
where

Center Manifold Reduction.

In order to reduce the order of the system, we apply the center manifold reduction [11–13]. We follow the steps defined in Ref. [13] for maps: Eq. (34) is already in Jordan form, so we can define the local center manifold in the form
(35)
where h is a polynomial function that satisfies Eq. (34) only for small values of ξ and η. The cubic terms are neglected, since they would appear only in higher order terms. We substitute Eq. (35) in the first and the second equation of Eq. (34), and then these two new equations and Eq. (35) into the third equation of Eq. (34). Collecting terms with the same power order we find
and solving the linear system we have
The system, in the vicinity of the bifurcation, consists of the first two equations of Eq. (34) where γ is substituted by the local center manifold. Substituting in Eq. (34) the coefficients phk as in Eq. (24), considering that at the bifurcation under study p=1/4 and d=5/8 and considering up to third order terms we have
(36)
(37)
Rewriting Eqs. (36) and (37) in matrix form we have
(38)

where ahk and bhk can be obtained from the nonlinear part of Eqs. (36) and (37).

Transformation to Normal Form.

In order to transform system (38) into its normal form in the vicinity of the bifurcation, we follow [13] and [14]: we rewrite the system in complex form
(39)
where

and ν,αijC.

The values of the coefficients are
(40)
(41)
We apply the nonlinear near identity transformation
(42)
to Eq. (39), where
It is possible to transform out the nonlinear terms by imposing
where α˜30, α˜12 and α˜03 are the coefficients of the cubic terms, modified because of the transformation of the second order terms, which is effecting them. In order to eliminate the term related to α21 we should have c21=-α˜21/(ν(1-νν¯)), that has no physical sense. So we let c21=0 and the normal form yields
(43)
(44)
where
(45)
The only nonlinear term left is related to the resonance μ1=μ12μ2 or μ2=μ22μ1 (at the bifurcation μ1μ2=1). Considering Eqs. (41) and (45), β21 can be found and it is equal to
(46)
Then we rewrite Eqs. (43) and (44) in the form
(47)
(48)
where
(49)

Bifurcation Diagram.

Considering that |v|2=vv¯ we transform Eqs. (47) and (48) into
(50)
where
(51)
Neglecting the sixth degree term, the solution for the radius of |v| of the invariant circle assumes the form
(52)
Considering Eq. (52), the type of bifurcation occurring depends strictly on the sign of δ, that is δ<0u02<3/5. The type of bifurcation does not depend on the sampling time τ, but only on the asymmetry of the control force. As already defined, u0=tanq0, so

moreover, according to the graph in Fig. 4, we must have |q0|<π/2 in order to have both positive and negative control force. Based on this result, the bifurcation under study is supercritical if the saturation of the force is symmetric (q0=0) or close to symmetry (|q0|<0.659) (Fig. 5), while it can become subcritical if the saturation is strongly asymmetric (|q0|>0.659).

In Fig. 5, the comparison between the analytical and the numerical bifurcation diagram is represented, as well. The matching around the critical value pcr=0.25 is very good. Due to the asymmetry, the periodic attractor is asymmetric itself; the bifurcation curves in Fig. 5 refer to the maximum absolute value of the mass displacement during its periodic motion.

From the analysis carried out in this section, if q0>0.659 the bifurcation will be subcritical, so for p<0.25 there will be a stable focus, i.e., the origin, and a branch of unstable solutions, while for p>0.25 there will be only an unstable focus that is the origin again. Since the system is bounded there must be always a stable attractor. In the case of subcritical bifurcation (q0>0.659) there should be at least another attractor that cannot be found with the local procedure carried out in this chapter since it has probably a large amplitude, so it is not caught by the accomplished approximation of the control force. In Sec. 6, the results of numerical simulations will provide a more detailed description of the system dynamics.

Numerical Simulation of the Full System

We carry out numerical simulations of the map corresponding to the system in Eq. (21). The aim of these simulations is to verify the results obtained from the analytical treatment and to find other attractors if they exist.

Most of the numerical analysis is made fixing the dimensionless differential gain d=5/8 and varying the parameters p and q0. Through continuation analysis and systematic numerical simulations, we obtain the response chart in Fig. 6, that shows the existing attractors of the map, given by Eqs. (21), (26) and (27), in the parameter space q0 and p. Attention is focused to a neighborhood of the maximal stability value of p and to a range of q0 where the Neimark-Sacker bifurcation is subcritical. Periodic and chaotic attractors are seen to occur both below p=0.25 and, mostly, above it. The region 5 can be extended down to q0=0. The dashed line limiting the chaotic region indicates that within it also non chaotic attractors may exist and the chaotic attractor may disappear in some parts. Regions 2 and 5 present more than one periodic attractor in some subregions, with up to three of them coexisting with each other.

Bifurcation Analysis.

In order to analyze the bifurcations outlined in Fig. 6, we plotted the bifurcation curves by varying one parameter per time, crossing all the boundaries in the diagram of Fig. 6.

The bifurcation occurring when crossing the line between regions 3 and 5 is represented by the diagrams in Fig. 7. According to the diagram in Fig. 6, the bifurcation under study is supercritical up to q00.65, which is in good accordance with the analytical results.

The typical bifurcation diagram ensuing from a subcritical Neimark-Sacker (Fig. 8) marks the passage from region 2 to region 5. It can be divided into three regions: for p<0.2492 we are in region 3 of Fig. 6, with the only attractor being the stable focus at the origin. For 0.2492<p<0.25 we are in region 2, with three coexisting solutions, the stable focus in the origin, a stable periodic attractor and a periodic repellor in between. For p>0.25 we are in region 5, where the latter disappears and the stable focus in the origin becomes unstable. This diagram has been obtained plotting for each value of p, the maximum displacement (in absolute value) after the transient has faded away, by sweeping the control parameter up and down in order to find possibly coexisting attractors.

It is clear that if we are in region 3, with p<0.249207 and we increase p, the mass will stay at the origin until we reach p=0.25, where it will start to oscillate and jump from point C to D; increasing the parameter p even more, the amplitude of vibration will progressively increase. If instead we start with a proportional gain p>0.25, the mass will start oscillating, with a progressive reduction of the amplitude of vibration as we decrease p until reaching point B at p=0.2492, then the mass will suddenly settle down onto the origin at point A. Furthermore, Fig. 8 shows a good matching between the numerical and the analytical results, around p = 0.25. While the amplitude increases, the matching becomes poorer because of the influence of the terms higher than the third.

The bifurcation occurring while crossing the border between regions 3 and 2, through a change in the parameter q0—which affects the asymmetry of the control force (see Fig. 10)—is represented in Fig. 9. For q0<0.668 we are in region 3 of the diagram in Fig. 6, as said before, where the only attractor is the stable focus in the origin. For q0>0.668 we are in region 2, with three coexisting solutions, as already explained. For higher values of q0 the motion becomes chaotic. The letters A and B of Fig. 9 correspond to those in Fig. 8. Also the bifurcation diagram in Fig. 9, shows an excellent agreement between the numerical and the analytical results, regarding the branch of unstable periodic solutions.

For high values of q0, the unstable branch tends to go close to the origin; thus reducing the robustness of the stable focus; at the same time the amplitude of the periodic attractor strongly increases. The combination of these two factors makes this region very dangerous from an engineering point of view, since, even if the parameters p and d guarantee stability, small perturbations to the equilibrium position can cause large oscillations.

Some parts of the stable branch are broken by other lines corresponding to a phase locking of the underlying continuous time system (see later) or to the coexistence of more than one attractor. Increasing the value of q0 the stable branch tends to be more and more rugged, highlighting a sequence of closer and closer bifurcations when approaching the chaotic region.

Transition to Chaos.

Figure 11 shows the variation of the attractors of the map given by Eqs. (21), (26) and (27) for different values of q0, i.e., for increasing asymmetry of the control force. In the figure, the displacement and the velocity are plotted for each time interval τ, so they can be considered as stroboscopic Poincaré maps of the underlying continuous time system. The figures on the top of Fig. 11 show a closed line, i.e., a periodic motion of the map, which corresponds to a quasi-periodic motion of the underlying continuous time system. For some specific parameter values, the continuous closed line of the Poincaré map is reduced to a finite number of points (q0=0.9), this corresponds to a phase locking of the different frequencies of the quasi-periodic motion. In the bifurcation diagrams of the map (Figs. 7–9), the phase locking shows up as an almost horizontal line breaking the continuous one. Figure 12 shows the regions, in the q0 and p parameter space, where there is a phase locking. These regions assume the structure of tongues, indeed similar structures are usually called Arnold tongues [15].

In Fig. 11, it is also visible how the attractor moves to the left and becomes more and more asymmetric for larger values of q0. Increasing the value of q0 even more, the motion is not periodic anymore and the attractor looks like being a strange attractor with a fractal structure, so there is a likely chaotic motion. The mechanism leading to chaos should be analyzed in more details, but this is out of the scope of this paper.

In order to verify if the attractor of the stroboscopic map (Eqs. (21), (26) and (27)), for the fixed values q0=0.97, p=0.2501Q0=10m/s2 and τ=0.01 s (Figs. 13(a) and 13(b)) is actually chaotic, we calculated a positive Lyapunov exponent [16,17], with the result being shown in Fig. 13(c). We also calculated its fractal dimension, using the Correlation Dimension [18,19], the result for this particular case being 1.4, and we plotted the frequency spectrum (Fig. 13(d)), which shows a broad-band response.

Conclusions and Suggestions for Further Research

In this paper, we first defined analytically the stability region of a single-DoF model robotic arm subjected to a digital PD control force, then we analyzed the bifurcation occurring at the stability border both analytically and numerically. The two results are perfectly matching in the case of a symmetric saturation of the control force, and both of them show the existence of a supercritical Neimark-Sacker bifurcation. In the case of an asymmetric saturation, above a certain, well-matched, level of asymmetry, they both show the existence of a subcritical Neimark-Sacker bifurcation, which may compromise the robustness of the stable solution in proximity of the stability border. Numerical simulations show the existence of chaotic motions for a high asymmetry of the saturation of the control force, established after a series of bifurcations. The onset and the characterization of chaotic motions deserve further investigations, but, from an engineering point of view, before the chaotic motion appears the periodic attractor is already too large to be accepted in a real application, so the analysis of the chaotic motion is not essential in this respect.

Improvements to the used model may give a better plausibility to the analysis; for example by considering also the dry friction acting on the mass and the finite stiffness of the arm of the robot.

References

1.
Gorinevsky
,
D. M.
,
Formalsky
,
A. M.
, and
Schneider
,
A. Y.
,
1997
,
Force Control of Robotics Systems
,
CRC Press LLC
,
Boca Raton, FL
.
2.
Siciliano
,
B.
, and
Villani
,
L.
,
1999
,
Robot Force Control
,
Kluwer
,
Dordrecht
.
3.
Natale
,
C.
,
2003
,
Interaction Control of Robot Manipulators
,
Springer
,
New York
.
4.
Colgate
,
E.
, and
Schenkel
,
G. G.
,
1997
, “
Passivity of a Class of Sampled-Data Systems: Application to Haptic Interfaces
,”
J. Rob. Syst.
14
(
1
), pp.
37
47
.
5.
Gil
,
J. J.
,
Avello
,
A.
,
Rubio
,
A.
, and
Florez
,
J.
,
2004
, “
Stability Analysis of a 1 DOF Haptic Interface Using the Routh-Hurwitz Criterion
,”
IEEE Trans. Control Syst. Technol.
12
, pp.
583
588
.10.1109/TCST.2004.825134
6.
Kovacs
,
L. L.
,
Kovecses
,
J.
, and
Stepan
,
G.
,
2008
, “
Analysis of Effects of Differential Gain on Dynamic Stability of Digital Force Control
,”
Int. J. Nonlinear Mech.
43
, pp.
514
520
.10.1016/j.ijnonlinmec.2008.04.002
7.
Stepan
,
G.
,
2001
, “
Vibrations of Machines Subjected to Digital Force Control
,”
Int. J. Solids Struct.
,
38
, pp.
2149
2159
.10.1016/S0020-7683(00)00158-X
8.
Hamzi
,
B.
,
Barbot
,
J. P.
,
Monaco
,
S.
, and
Normand-Cyrot
,
D.
,
2001
, “
Nonlinear Discrete-Time Control of Systems with a Naimark-Sacker Bifurcation
,”
Syst. Control Lett.
,
44
(
4
), pp.
245
258
.10.1016/S0167-6911(01)00136-0
9.
Yaghoobi
,
H.
, and
Abed
,
E. H.
,
2003
, “
Local Feedback Control of the Naimark-Sacker Bifurcation
,”
Int. J. Bifurcation Chaos Appl. Sci. Eng.
,
13
(
4
), pp.
879
893
.10.1142/S0218127403006972
10.
Yuan
,
Z.
,
Hu
,
D.
, and
Huang
,
L.
,
2005
, “
Stability and Bifurcation Analysis on a Discrete-Time Neural Network
,”
J. Comput. Math.
,
177
, pp.
89
100
.10.1016/j.cam.2004.09.010
11.
Guckenheimer
,
J.
, and
Holmes
,
P.
,
1986
,
Nonlinear Oscillations, Dynamical Systems, and Bifurcation of Vector Fields
,
Springer
,
New York
.
12.
Nayfeh
,
A. H.
, and
Balachandran
,
B.
,
1995
,
Applied Nonlinear Dynamics
,
Wiley
,
New York
.
13.
Troger
,
H.
, and
Steindl
,
A.
,
1991
,
Nonlinear Stability and Bifurcation Theory
,
Springer
,
New York
.
14.
Szalai
,
R.
,
Stepan
,
G.
, and
Hogan
,
S. J.
,
2004
, “
Global Dynamics of Low Immersion High-Speed Milling
,”
Chaos.
14
(
4
), pp.
1069
1077
.
15.
Kuznetsov
,
Y. A.
,
1998
,
Elements of Applied Bifurcation Theory
,
Springer
,
New York
.
16.
Moon
,
F. C.
,
1992
,
Chaotic and Fractal Dynamics: an Introduction for Applied Scientists and Engineers
,
Wiley
,
New York
.
17.
Farmer
,
J. D.
,
Ott
,
E.
, and
York
,
J. A.
,
1983
, “
The Dimension of Chaotic Attractors
,”
Physica D
,
7
,
153
170
.10.1016/0167-2789(83)90125-2
18.
Grassberger
,
P.
, and
Procaccia
,
I.
,
1983
, “
Characterization of Strange Attractors
,”
Phys. Rev. Lett.
,
50
(
5
), pp.
346
349
.10.1103/PhysRevLett.50.346
19.
Hilborn
,
R. C.
,
2001
,
Chaos and Nonlinear Dynamics: an Introduction for Scientists and Engineers
,
Oxford University
,
New York
.