Stability analysis of admittance control
using asymmetric stiffness matrix

In-demand automation tasks often involve contact-rich tasks with varying mechanical constraints between the robot and the environment, making control program implementation complex. Studies on contact-rich tasks using the peg-in-hole benchmark have expanded automation technology, particularly for precision assembly. Compliance mechanisms, such as active and passive compliance, help robots to perform assembly tasks despite positional errors. Impedance and admittance control are common methods for active compliance [1, 2]. Contact-model-based [3, 4] and contact-model-free [5, 6] approaches are used for compliance control, with the latter often using machine learning methods such as learning from demonstration and learning from the environment with reinforcement learning.

Variable impedance control (VIC), a control strategy that adjusts a robot’s mechanical impedance for safe, versatile, and efficient interaction, has been actively studied [7]. VIC has shown promising results in various applications, such as manipulation, locomotion, haptic feedback, and rehabilitation. However, it faces several challenges, such as the design of robust and efficient actuators and sensors, as well as the development of effective control strategies for complex and dynamic tasks [8].

Methods that use passive compliance obtain the desired stiffness of the hand by using a device that combines passive mechanical elements such as springs. The remote center compliance (RCC) device is one of the most popular examples, which utilizes the inter-axis interference to induce motion against external force [9]. Instead of VIC, variable RCC (VRCC), which can set stiffness in a variable manner, has been studied to deal with complicated tasks [10, 11].

RCC has the advantage of not being limited by control bandwidth because it is composed of mechanical elements. On the other hand, when compliance is reproduced by control, even ranges that cannot be reproduced mechanically can be artificially reproduced with a single parameter setting. In other words, the compliance control has a wider range over which motion can be induced. However, how to set the parameters under an assurance of stability is an important issue.

A passivity-based approach is the most popular candidate to ensure the stability of VIC [12]. It has also been extended to the time-varying admittance controller [13]. Kronander and Billard proposed stability conditions for the VIC scheme that works for different levels of stiffness and damping, using a modified Lyapunov function [14]. Later, Sun et al. expanded on this idea and proposed new constraints that ensure both exponential stability and boundedness of the robot’s position, velocity, and acceleration [15]. Spyrakos-Papastavridis et al. proposed a Passivity-Preservation Control (PPC) that enables the implementation of stable VIC and provided joint and Cartesian space versions of the PPC controller permit the intuitive definition of interaction tasks [16]. In the study of parameter design for variable admittance control, Kim and Yang performed a stability analysis based on the root locus [17]. The root locus has also been used in the parameter design for drill admittance control [18].

While many stability analyses of compliance control have been conducted, only passivity-based analysis has been pursued in multi-degree-of-freedom systems because the problem becomes complex due to the consideration of interference in other axes [19]. Previous studies on compliance control, such as [20] and [21], have indicated that introducing an asymmetric stiffness matrix can increase the robot’s capability and extend the range of motion induction. However, the main challenge is that the asymmetric part fails to satisfy passivity, making it impossible to ensure stability. In this study, we derived parameters for the stability limit of the asymmetric stiffness matrix based on the eigenvalues of the system matrix. We then demonstrate that stability can be guaranteed by setting the parameters using the root locus method.

The contributions of this paper can be summarized as follows.

1. 1.

   Spiral oscillations can be excited when the stiffness matrix is asymmetric. We have clarified the condition for such excitation. We also showed that the condition for the absence of spiral oscillation is equivalent to the condition that the eigenvalues of the stiffness matrix are real numbers.

2. 2.

   We have demonstrated that it is difficult to satisfy the passivity when the above conditions are not satisfied and therefore identified the conditions under which vibration converges using the root locus method.

This study discusses the stability of compliance control in a robot arm that deals with an asymmetric stiffness matrix. It is limited to admittance control and discusses the convergence stability of a planar admittance model, assuming the inner position control loop is stable. The feedback control system is shown in Fig. 1. It is an admittance control system with a minor feedback loop for position control. The type of position control is not specified, but a control system using a disturbance observer (DOB) is adopted in the experimental setup. The admittance model is based on the following equation:

$\displaystyle\mbox{\boldmath$F$}_{ext}=\mbox{\boldmath$M$}\ddot{\mbox{\boldmath$x$}}+\mbox{\boldmath$D$}\dot{\mbox{\boldmath$x$}}+\mbox{\boldmath$K$}\mbox{\boldmath$x$}$

(1)

This admittance model takes force $\mbox{\boldmath$F$}_{ext}$ as input and outputs the position $x$ and velocity $\dot{\mbox{\boldmath$x$}}$ obtained through numerical integration based on the acceleration $\ddot{\mbox{\boldmath$x$}}$ obtained from equation (1). While symmetric matrices are typically used for the stiffness matrix $K$ and damping matrix $D$ of the admittance model, this study considers the use of an asymmetric matrix for the stiffness matrix $K$. For simplicity, the damping matrix $D$ is assumed to be a diagonal matrix.

When an external force occurs, the direction in which the robot end effector is induced, and the convergence point are determined only by the stiffness matrix $K$ of the admittance model. Therefore, it is desirable to expand the selection range of the stiffness matrix $K$ when designing the relationship between external forces and robot operation. On the other hand, since the damping matrix $D$ does not affect the convergence point of the end effector, the priority of expanding the selection range is not high. Additionally, applying viscous force only in the direction in which velocity is generated increases the dissipation of energy and is reasonable for enhancing stability. Therefore, it is necessary to consider how to give the diagonal elements of the damping matrix $D$ when the non-diagonal elements of the stiffness matrix $K$ are given arbitrarily.

When the stiffness matrix $K$ is asymmetric, it is known that the effect of the stiffness matrix can be expressed as a superposition of the symmetric part $\mbox{\boldmath$K$}_{s}$, which creates a force field that converges at the equilibrium point [1], and the asymmetric part $\mbox{\boldmath$K$}_{a}$, which creates a force that rotates around the equilibrium point. An example of this is shown in Fig. 2. If the stiffness control is reproduced by a virtual spring, the symmetric part can be treated as generating an elliptical virtual potential energy field generated by strain. Similarly, virtual kinetic energy can be derived from the mass and velocity of the admittance model. Therefore, in the symmetric matrix, a control system design based on passivity can be achieved by deriving an energy function. However, when this is extended to an asymmetric matrix, the curl field generated by the asymmetric part is not a conservative force, and it may excite oscillations. Oscillations excited by the curl field are referred to as spiral oscillations. Since energy may increase due to spiral oscillation, passivity may not be satisfied.

This is demonstrated by the energy function. It is obtained as the integral of the product of velocity and force, which is power.

$\displaystyle V\!\!$

$\displaystyle=$

$\displaystyle\int\dot{\mbox{\boldmath$x$}}^{T}\mbox{\boldmath$F$}_{ext}dt$

(2)

By substituting $\mbox{\boldmath$K$}=\mbox{\boldmath$K$}_{s}+\mbox{\boldmath$K$}_{a}$ and (1) into (2),

	$\displaystyle V\!\!$	$\displaystyle=$	$\displaystyle\int\dot{\mbox{\boldmath$x$}}^{T}\mbox{\boldmath$M$}\ddot{\mbox{\boldmath$x$}}+\dot{\mbox{\boldmath$x$}}^{T}\mbox{\boldmath$D$}\dot{\mbox{\boldmath$x$}}+\dot{\mbox{\boldmath$x$}}^{T}\mbox{\boldmath$K$}_{s}\mbox{\boldmath$x$}+\dot{\mbox{\boldmath$x$}}^{T}\mbox{\boldmath$K$}_{a}\mbox{\boldmath$x$}$		(3)
		$\displaystyle=$	$\displaystyle\!\!\!\!\frac{1}{2}\dot{\mbox{\boldmath$x$}}^{T}\!\!\mbox{\boldmath$M$}\dot{\mbox{\boldmath$x$}}+\frac{1}{2}\mbox{\boldmath$x$}^{T}\!\!\mbox{\boldmath$K$}_{s}\mbox{\boldmath$x$}+\!\!\int\!\!\dot{\mbox{\boldmath$x$}}^{T}\!\!\mbox{\boldmath$D$}\dot{\mbox{\boldmath$x$}}dt+\!\!\int\!\!\dot{\mbox{\boldmath$x$}}^{T}\!\!\mbox{\boldmath$K$}_{a}\mbox{\boldmath$x$}dt$		(4)

Eq. (4) shows that the virtual energy of the admittance model is defined as a function divided into four terms: kinetic energy, stiffness energy, energy loss due to friction, and energy imparted by the asymmetric part. The first and second terms represent kinetic and potential energy, respectively. The third term represents energy dissipated by viscosity, while the fourth term represents energy imposed by the curl field. Since the first and second terms are the energy of conservative forces, if the third term is smaller than the fourth term, the energy function does not increase. However, the above equation depends on velocity, and when the velocity decreases, the dissipative term becomes smaller. Therefore, if there are periods during vibration where the velocity becomes small, passivity may not hold.

The fundamental principle of admittance control described in the previous section shows that it is generally difficult to establish passivity when the stiffness matrix is made asymmetric. Therefore, the stability will be analyzed based on the concept of root locus below. First, the admittance parameters are defined in the x-y plane as follows:

$$\mbox{\boldmath$K$}=\mbox{\boldmath$K$}_{s}+\mbox{\boldmath$K$}_{a}=\begin{bmatrix}k_{x}&k_{s}+k_{a}\\ k_{s}-k_{a}&k_{y}\end{bmatrix},\mbox{\boldmath$D$}=\begin{bmatrix}d&0\\ 0&d\end{bmatrix}$$

(5)

we define the state vector $x$ as $[x_{x},x_{y},\dot{x}_{x},\dot{x}_{y}]$, then its state equation can be obtained as follows:

	$\displaystyle\dot{\mbox{\boldmath$x$}}$	$\displaystyle=$	$A$$x$
	$\displaystyle\begin{bmatrix}\ddot{x}_{x}\\ \ddot{x}_{y}\\ \dot{x}_{x}\\ \dot{x}_{y}\end{bmatrix}$	$\displaystyle=$	$\displaystyle\begin{bmatrix}-\frac{d}{m}&0&-\frac{k_{x}}{m}&\frac{-k_{s}-k_{a}}{m}\\ 0&-\frac{d}{m}&\frac{-k_{s}+k_{a}}{m}&-\frac{k_{y}}{m}\\ 1&0&0&0\\ 0&1&0&0\end{bmatrix}\begin{bmatrix}\dot{x}_{x}\\ \dot{x}_{y}\\ x_{x}\\ x_{y}\end{bmatrix}$		(6)

Here, $k_{x},k_{y},k_{s},k_{a},d,m\in\mbox{\boldmath$R$}$. $d$ and $m$ are the damper and mass of the admittance model. $k_{x}$ and $k_{y}$ represent the diagonal elements of the stiffness matrix in the x and y directions, respectively, and $k_{s}$ and $k_{a}$ represent the symmetric and asymmetric components of the non-diagonal terms, respectively. Moreover, we define parameters as follows: $d_{m}=\frac{d}{m},k_{xm}=\frac{k_{x}}{m},k_{ym}=\frac{k_{y}}{m},k_{sm}=\frac{k_{s}}{m},k_{am}=\frac{k_{a}}{m}$, then we can write $A$, the system matrix of the state equation, as follows:

$A$

$\displaystyle=$

$\displaystyle\begin{bmatrix}-d_{m}&0&-k_{xm}&-k_{sm}-k_{am}\\ 0&-d_{m}&-k_{sm}+k_{am}&-k_{ym}\\ 1&0&0&0\\ 0&1&0&0\end{bmatrix}$

(7)

The eigenvalues of the system matrix $A$ are obtained as follows:

$$\!\!\lambda\!=\!\!\frac{\!\!-d_{m}\!\!\pm\!\!\sqrt{\!d_{m}^{2}\!\!-\!\!2(k_{xm}\!\!+\!\!k_{ym}\!\!\pm\!\!\sqrt{\!\!-4k_{am}^{2}\!\!+\!\!4k_{sm}^{2}\!\!+\!\!(k_{xm}\!\!-\!\!k_{ym})^{2}\!})}}{2}$$

(8)

The poles of the transfer function of this admittance model have the same values as the eigenvalues of the system matrix $A$. Therefore, the root locus can be obtained by sweeping a parameter of matrix $A$ from 0 to infinity. When the non-diagonal element $k_{am}$ of the stiffness matrix is zero, i.e., when the stiffness matrix is symmetric, the root locus for varying the damping ratio $d_{m}$ from 0 to infinity is shown in Fig. 3. For a symmetric matrix, the roots of the natural vibrations, which depend on the stiffness matrix, appear on the imaginary axis, and as the viscous damping ratio $d_{m}$ increases, the roots move in the negative real axis direction. At the same time, the imaginary part of the root decreases. The damping ratio $d_{m}$ where the imaginary part of the two conjugate eigenvalues shrinks to zero is called the critical damping ratio. The contents of the square root of (8) are then zero, hence the critical damping ratio is derived by the following equation.

$\displaystyle d_{m}\!\!=\!\sqrt{2\left(k_{xm}\!+\!k_{ym}\!\!\pm\!\!\sqrt{\!-4k_{am}^{2}\!\!+\!4k_{sm}^{2}\!\!+\!(k_{xm}\!\!-\!k_{ym})^{2}}\right)}$

(9)

When the value is equal to the critical damping ratio, the roots coincide on the real axis. When the viscous damping ratio $d_{m}$ increases beyond this value, the roots separate and move to the left and right along the real axis.

Fig. 3(a) shows the result when $k_{xm}=k_{ym}=100$ and $k_{sm}=0$, and the stiffness matrix is diagonal. On the other hand, Fig. 3(b) shows the result when $k_{xm}=k_{ym}=100$ and $k_{sm}=40$. In other words, the stiffness matrix is not diagonal and its eigenvalues are 140 and 60. Fig. 3(c) shows the result for the case of a diagonal stiffness matrix with the same eigenvalues as in Fig. 3(b), with $k_{xm}=140,k_{ym}=60,k_{sm}=0$. The initial position of the root locus changes depending on the eigenvalues of the stiffness matrix, as shown in the comparison of Figs. 3(a), (b), and (c).

In Fig. 3(a), since the two eigenvalues coincide, the root locus starts from the points $\pm\sqrt{k_{xm}}j$ on the imaginary axis, while in Fig. 3(b), the root locus is split into two upper and lower branches due to the different eigenvalues. Figs. 3(b) and (c) also show that the root locus coincides when the eigenvalues are the same. Eq. (1) deals with a generalized stiffness matrix in two dimensions, but to simplify the derivation of the equation, $k_{sm}$ is set to 0, and stability analysis is performed below. For the case where $k_{sm}$ is not 0, the same root locus can be obtained by deriving its eigenvalues and substituting them into the diagonal elements.

The root locus of the asymmetric matrix is shown in Fig. 4, where four roots appear symmetrically across the imaginary axis when $d_{m}$ is 0. As $d_{m}$ is gradually increased from 0, the four roots move negatively along the real axis while the imaginary axis component decreases. As $d_{m}$ is further increased, two of the four roots change their orientation toward the positive direction of the real axis, and the other two roots move in opposite directions, one positive and one negative. In this case, the imaginary axis component remains slightly. This means that in the case of an asymmetric matrix, oscillations are excited due to the presence of the curl field. However, there is a condition for the parameter that excites the oscillation. To ensure that all eigenvalues $\lambda$ are real, it is necessary for the contents of the square root in the second term of the numerator in (8) to be positive, that is,

$$d_{m}^{2}\!-2\!\left(\!k_{xm}\!+\!k_{ym}\!\pm\!\!\sqrt{\!-4k_{am}^{2}\!+4k_{sm}^{2}\!+\!\left(k_{xm}\!-\!k_{ym}\right)^{2}}\right)\!\!>\!0.$$

(10)

And it is necessary to have a real value of $d_{m}$ in (10). To ensure this, the square root term in (10) must be positive. The condition is expressed by the following inequality.

$$-4k_{am}^{2}+4k_{sm}^{2}+(k_{xm}-k_{ym})^{2}>0$$

(11)

Consequently, condition for the asymmetric element $k_{am}$ is determined as follows:

$$|k_{am}|<\sqrt{k_{sm}^{2}+1/4(k_{xm}-k_{ym})^{2}}$$

(12)

If the above inequality is satisfied, the oscillation will no longer occur when $d_{m}$ is greater than the critical damping. The eigenvalues of the stiffness matrix $K$ are derived by the following equation:

$$\lambda_{K}=\frac{1}{2}\left(k_{xm}+k_{ym}\pm\sqrt{-4k_{am}^{2}+4k_{sm}^{2}+\left(k_{xm}-k_{ym}\right)^{2}}\right)$$

(13)

The condition for the stiffness matrix to obtain critical damping is equivalent to the matrix having real eigenvalues, as shown by the fact that the condition for (13) to have real solutions is the same as that for (12).

Next, we will derive the conditions for the admittance parameters that ensure stability. If the viscous damping coefficient $d_{m}$ is increased to a certain value or more, all the roots will move toward the left half of the complex plane, making it possible to stabilize the control system. From (8), in the case of an asymmetric matrix,

$\displaystyle k_{xm}+k_{ym}\pm\sqrt{(-4k_{am}^{2}+(k_{xm}-k_{ym})^{2}}<0$

(14)

The above condition results in the occurrence of positive eigenvalues and unstable poles. However, if $-4k_{am}^{2}+(k_{xm}-k_{ym})^{2}>0$, there are no solutions that satisfy (14) for $k_{xm}>0,k_{ym}>0$. Therefore, the condition for unstable poles to occur within the range of $-4k_{am}^{2}+(k_{xm}-k_{ym})^{2}\leq 0$ needs to be investigated. To simplify the (8), $\alpha,\beta\in\mbox{\boldmath$R$}$ are defined as follows:

$\displaystyle\alpha=\!d_{m}^{2}\!-\!2(k_{xm}\!+\!k_{ym}),\beta=\!\!\sqrt{16k_{am}^{2}\!-\!4(k_{xm}\!-\!k_{ym})^{2}}.$

(15)

Then, the eigenvalues are expressed as follows:

$\displaystyle\lambda=\frac{1}{2}(-d_{m}\pm\sqrt{\alpha\pm\beta j})$

(16)

Let the real part and imaginary part of $\sqrt{\alpha\pm\beta j}$ be $a$ and $b$, respectively, then

$\displaystyle(a+bj)(a+bj)=a^{2}-b^{2}+2abj=\alpha+\beta j$

(17)

Consequently, $\alpha=a^{2}-b^{2}$, and $\beta=2ab$ are obtained. By the following equation: $\alpha+\sqrt{\alpha^{2}+\beta^{2}}=2a^{2}$, real part $a$ is expressed as below:

$\displaystyle a=\pm\sqrt{\frac{\alpha+\sqrt{\alpha^{2}+\beta^{2}}}{2}}$

(18)

In order to avoid unstable poles, it is necessary for all eigenvalues to be negative, so we only consider the case where $a=\sqrt{\frac{\alpha+\sqrt{\alpha^{2}+\beta^{2}}}{2}}$.

$\displaystyle-d_{m}+\sqrt{\frac{\alpha+\sqrt{\alpha^{2}+\beta^{2}}}{2}}<0$

(19)

The condition that satisfies the (10) is derived below. First, consider the typical case where the equation is simplified, which is $k_{xm}=k_{ym}$. Substituting $\alpha=d_{m}^{2}-4k_{xm}$, $\beta=4k_{am}$ into (19) and the following equation is obtained:

$\displaystyle d_{m}>\sqrt{\frac{k_{am}^{2}}{k_{xm}}}$

(20)

If the value of $k_{am}$ is greater than $k_{xm}$, then the lower limit of $d_{m}$ also increases. These parameters are determined by the mass ratio, and to convert them into actual admittance parameters, we substitute $k_{xm}=k_{x}/m,k_{am}=k_{a}/m$, and $d_{m}=d/m$ into the above equation, yielding the following equation:

$\displaystyle d>\sqrt{\frac{k_{a}^{2}m}{k_{x}}}$

(21)

The next is to consider the case where $k_{ym}>k_{xm}>0$ to further generalize (20) and (21). The condition is still given by (15), but now we substitute the more generalized (10). It is difficult to algebraically derive the condition for all eigenvalues to be negative. However, if (20) is satisfied, then all eigenvalues will be negative.

The proof is shown below. If (19) is satisfied and without changing the values of the other parameters, let us increase the value of $k_{ym}$ from $k_{xm}$. If $a$ in (18) does not increase during this process, (15) will always be satisfied. When differentiating (18) with respect to $\alpha$ and $\beta$, we obtain the following expressions:

	$\displaystyle\frac{\partial a}{\partial\alpha}$	$\displaystyle=$	$\displaystyle\frac{1+\frac{a}{\sqrt{\alpha^{2}+\beta^{2}}}}{2\sqrt{2}\sqrt{\alpha+\sqrt{\alpha^{2}+\beta^{2}}}}$		(22)
	$\displaystyle\frac{\partial a}{\partial\beta}$	$\displaystyle=$	$\displaystyle\frac{\beta}{2\sqrt{2}\sqrt{\alpha^{2}+\beta^{2}}\sqrt{\alpha+\sqrt{\alpha^{2}+\beta^{2}}}}$		(23)

when $\beta>0$ and $k_{ym}>k_{xm}$, both partial derivatives are positive. As (15) shows, increasing $k_{ym}$ always results in a decrease in both $\alpha$ and $\beta$, so the value of $a$ in (18) must always decrease. Therefore, if (20) is satisfied, then for all $k_{ym}>k_{xm}>0$, all eigenvalues are always negative. Replacing $k_{xm}$ and $k_{ym}$ in (20) as follows:

$\displaystyle d_{m}>\sqrt{\frac{k_{am}^{2}}{k_{ym}}}$

(24)

By performing the same equation expansion as above, stability can also be demonstrated for the case where $k_{xm}>k_{ym}>0$.

In this experiment, the robot’s performance was evaluated based on the admittance parameter. The robot was made to move vertically towards a metallic board and follow the surface in the y-axis using the displacement caused by the admittance control. Symmetric and asymmetric stiffness matrices were chosen, and damping parameters were set to $d=18$ which was higher than the threshold ($d=17.889$). The results shown in Fig. 9 revealed that after contacting the board at $t=27$ seconds, the robot moved towards the negative direction in the y-axis. There was no noticeable movement in the x-axis due to the design of the non-diagonal element of the stiffness matrix. In the case of the symmetric matrix, the difference between the desired displacement $y_{d}$ and the actual displacement $y$ was 8 mm. Due to the setting of $k_{zy}=400$, there was slight displacement in the z-axis, which led to a lower contact force $f_{z}$. However, in the case of the asymmetric (triangular) matrix, the difference between the desired displacement $y_{d}$ and the actual displacement $y$ was 6 mm, which is shorter than that of the symmetric matrix. Additionally, the contact force $f_{z}$ was maintained higher because the displacement mainly occurred in the y-axis.

This paper discusses the stability of admittance control of a robot arm handling an asymmetric stiffness matrix. Although the introduction of asymmetric elements does not guarantee passivity, it is possible to derive a range of parameters for which stability is guaranteed from the eigenvalues of the system matrix. In addition, the introduction of asymmetric elements may excite helical oscillations, so that critical damping may not be obtained for any given damping parameter, depending on the setting of the stiffness matrix. In this study, we focused on the fact that the spiral oscillation is caused by the imaginary component of the system matrix and clarified the condition of the asymmetric matrix for which critical damping can be obtained. From the simulation and experimental results, it is confirmed that the asymmetric matrix setting method is effective to obtain critical damping. The analysis was performed on a 2-dimensional plane, but future work should extend the theory to a wider range of conditions. The theory can be generalized to systems with three or more degrees of freedom and by utilizing the concept of energy tanks, conduct stability analysis similar to a passive system. The theory may also be extended to variable stiffness control.