## 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\xb7$ of the mass, and elaborates a current signal that is sent to the DC motor.

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.

which is piecewise constant.

## Stability Analysis

### Discrete Mapping.

*j*+ 1, is equal to 1. So in matrix form we have

### 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 $|\mu 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 same stability condition can be derived using the Laplace transformation $L$, which is illustrated in Fig. 2.

the $kth$ column of the transformation matrix contains the binomial coefficients of $(\sigma +1)4-k\xd7(\sigma -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.

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\tau 2/m=1/4$ and $D\tau /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 $\mu 1,2=e\xb1i\alpha $. According to Floquet theory in this circumstance there is a Neimark-Sacker bifurcation [11].

### Vibration Frequencies.

## Nonlinear Analysis

### Nonlinear System.

*m*. Eq. (1) becomes

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).

### Jordan Normal Form.

*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 $p\u2243pcr=1/4$ are

### Center Manifold Reduction.

*h*is a polynomial function that satisfies Eq. (34) only for small values of $\xi $ and $\eta $. 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

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

### Transformation to Normal Form.

and $\nu ,\alpha ij\u2208C$.

### Bifurcation Diagram.

moreover, according to the graph in Fig. 4, we must have $|q0|<\pi /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 $q0\u22480.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 $\tau $, 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.2501$$Q0=10\u2003m/s2$ and $\tau =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.