Robot control based on motor primitives: A comparison of two approaches

Abstract

Motor primitives are fundamental building blocks of a controller which enable dynamic robot behavior with minimal high-level intervention. By treating motor primitives as basic “modules,” different modules can be sequenced or superimposed to generate a rich repertoire of motor behavior. In robotics, two distinct approaches have been proposed: Dynamic Movement Primitives (DMPs) and Elementary Dynamic Actions (EDAs). While both approaches instantiate similar ideas, significant differences also exist. This paper attempts to clarify the distinction and provide a unifying view by delineating the similarities and differences between DMPs and EDAs. We provide nine robot control examples, including sequencing or superimposing movements, managing kinematic redundancy and singularity, control of both position and orientation of the robot’s end-effector, obstacle avoidance, and managing physical interaction. We show that the two approaches clearly diverge in their implementation. We also provide a real-robot demonstration to show how DMPs and EDAs can be combined to get the best of both approaches. With this detailed comparison, we enable researchers to make informed decisions to select the most suitable approach for specific robot tasks and applications.

1. Introduction

One of the major challenges of robotics is to generate complex motor behavior that can match that of humans (Billard and Kragic, 2019). Numerous approaches have been developed to address this problem, including applied nonlinear control (Chung and Slotine, 2009; Slotine and Li, 1991), optimization-based approaches (Kuindersma et al., 2016; Posa et al., 2014), and machine learning algorithms (Haarnoja et al., 2018; Lillicrap et al., 2015; Schulman et al., 2015).

Among these methods, several advances have been made based on “motor primitives” (Hogan and Sternad, 2012; Ijspeert et al., 2013; Schaal, 1999). The fundamental concept originates from human motor control research, where complex motor behavior of biological systems appears to be generated by a combination of fundamental building blocks known as motor primitives (Bizzi and Mussa-Ivaldi, 2017; Flash and Hochner, 2005; Hogan and Sternad, 2012; Mussa-Ivaldi and Bizzi, 2000; Wolpert et al., 2001).

The concept of control based on motor primitives dates back at least a century (Sherrington, 1906), with a number of subsequent experiments providing support for its existence in biological systems. Sherrington was one of the first to suggest “reflex” as a fundamental element of complex motor behavior (Elliott et al., 2001; Sherrington, 1906). Sherrington proposed that reflexes can be treated as basic units of motor behavior that when chained together produce more complex movements (Clower, 1998). First formalized by Bernstein (Bernstein, 1935; Latash, 2021), “synergies” have also been suggested as a motor primitive to account for the simultaneous motion of multiple joints or activation of multiple muscles (Giszter et al., 1993; d’Avella et al., 2003; Giszter and Hart, 2013).

Along with reflex and synergy, kinematic motor primitives (Giszter, 2015) such as “rhythmic” and “discrete” movements have also been suggested as motor primitives (Hogan and Sternad, 2007, 2012; Huh and Sejnowski, 2015; Park et al., 2017; Schaal et al., 2004; Viviani and Flash, 1995). Rhythmic movements (e.g., locomotion) are phylogenetically old motor behaviors found in most biological species (Ronsse et al., 2009; Schaal et al., 2004). Central Pattern Generators (Brown, 1911, 1912; Marder and Bucher, 2001), which are specialized neural circuits for generating rhythmic motor patterns, have been identified in biological systems (Hogan and Sternad, 2012). In comparison, discrete movements (e.g., goal-directed reaching movements) are phylogenetically young motor behaviors, particularly observed in primates with developed upper extremities (Ronsse et al., 2009). Neural imaging studies have effectively ruled out the hypothesis that rhythmic (arm) movements result from the concatenation of discrete movements (Schaal et al., 2004), further supporting that rhythmic and discrete movements constitute distinct classes of primitives.

Recently, there is growing evidence that neural circuits responsible for maintaining posture are distinct from those that control movement. Hence, “stable postures” may be considered to be a distinct class of motor primitives (Jayasinghe et al., 2022; Shadmehr, 2017).

Motor primitives have also been applied to robotics. Two distinct approaches exist: Dynamic Movement Primitives (DMPs) (Schaal, 2006; Ijspeert et al., 2013; Saveriano et al., 2023) and Elementary Dynamic Actions (EDAs) (Hogan, 2017; Hogan and Sternad, 2012, 2013). ¹ The key idea of these approaches is to formulate motor primitives as “attractors” (Hogan and Sternad, 2012; Ijspeert et al., 2013; Saveriano et al., 2023). An attractor is a prominent feature of nonlinear dynamical systems, defined as a set of states toward which the system tends to evolve. Its type ranges from relatively simple ones such as (stable) “point attractors” and (stable) “limit cycles,” to “strange attractors” (Strogatz, 2018) such as the “Lorenz attractor” (Lorenz, 1963), “Rössler attractor” (Rössler, 1976), and others (Sprott, 2014; Tam et al., 2008).

One of the key benefits of using motor primitives is that it enables highly dynamic behavior of the robot with minimal high-level intervention (Hogan and Sternad, 2012). As a result, the complexity of the control problem can be significantly reduced. For instance, by formulating discrete (respectively rhythmic) movement as a stable point attractor (respectively limit cycle), problem of generating the movement reduces to learning the parameters of the corresponding attractor (Schaal et al., 2007). Another important consequence is that it provides a modular structure of the controller. By treating motor primitives as basic “modules,” learning motor skills happens at the level of modules which provides adaptability and flexibility for robot control.

Since DMPs and EDAs stem from the theory of motor primitives, both approaches share the same philosophy. Nevertheless, significant differences exist such that their implementations diverge. In the opinion of the authors, this has not yet been sufficiently emphasized. An in-depth review that elucidates the similarities and differences between the two approaches may be beneficial to the robotics community. Furthermore, a method that integrates the strength of both motor-primitives approaches could be beneficial to solve a wide range of control tasks.

In this paper, we provide a comprehensive review of motor primitives in robotics, focusing specifically on the two distinct approaches—DMPs and EDAs. We first delineate the similarities and differences of both approaches by presenting nine extensively used robotic control examples (Section 3). ² We show that (Table 1):

• Both approaches use motor primitives as basic building blocks to parameterize the controller (Section 2). DMPs consist of a canonical system, nonlinear forcing terms, and transformation systems. EDAs consist of submovements, oscillations, and mechanical impedances.

• For torque-controlled robots, DMPs require an inverse dynamic model of the robot, whereas EDAs do not impose this requirement (Section 3.1).

• With an inverse dynamics model, DMPs can achieve perfect tracking, both in task-space and joint-space (Section 3.2, 3.3). Imitation Learning enables DMPs to learn and track trajectories of arbitrary complexity (Section 2.1.4, 5.2.2). Online trajectory modulation of DMPs enables achieving additional control objectives such as obstacle avoidance (Section 3.5), thereby providing advantages over spline methods. For tracking control with EDAs, an additional method for calculating an appropriate virtual trajectory and mechanical impedance to which it is connected (Section 2.2.4) is required (Section 2.2, 3.2, 3.3).

• To control the position and/or orientation of the robot’s end-effector, DMPs require additional control methods to manage kinematic singularity and kinematic redundancy (Section 3.3, 3.9, 3.10). In contrast, for EDAs, stability near (and even at) kinematic singularity can be ensured (Section 3.3). Kinematic redundancy can be managed without solving the inverse kinematics (Section 3.9, 3.10).

• Both approaches provide a modular framework for robot control. However, the extent of modularity and its practical implications differ between the two approaches. A clear distinction appears when combining multiple movements. For DMPs, discrete and rhythmic movements are represented by different DMPs (Section 2.1.1, 2.1.2). Hence, the two different DMPs cannot be directly superimposed to generate a combination of discrete and rhythmic movements (Section 3.7). Multiple discrete movements are generated by modifying the goal position of the previous movement (Section 3.8). While the weights of the nonlinear forcing terms learned from Imitation Learning (Section 2.1.4) can be reused, the weights of different DMPs cannot be simply combined.

For EDAs, modularity exists both at the level of kinematics and mechanical impedances (Section 2.2.4). For the former, multiple movements can be directly superimposed at the kinematic level, which provides a greater degree of simplicity to generate a rich repertoire of movements, for example, combining discrete and rhythmic movements (Section 3.7), or sequencing discrete movements (Section 3.8). For the latter, the learned mechanical impedances can be reused and combined by simple superposition. Superposition of mechanical impedances enables a “divide-and-conquer” strategy, where complex tasks can be broken down into simpler sub-tasks. This modular property of mechanical impedances simplifies multiple control tasks, for example, obstacle avoidance (Section 3.5) and managing kinematic redundancy (Section 3.9, 3.10).

Table 1.

Summary of the main differences between dynamic movement primitives (DMPs) and elementary dynamic actions (EDAs).

	Dynamic movement primitives (DMPs)	Elementary dynamic actions (EDAs)
Controller parameterization	Canonical System, Nonlinear Forcing Terms, and Transformation Systems	Submovements, Oscillations, and Mechanical Impedances
Robot compatibility	Position- and torque-commanded* robots	Torque-commanded robots only
Tracking control	High tracking performance with an accurate robot model	Non-negligible tracking error without additional optimization, but without the need for an accurate robot model
	Direct command for position control but requires inverse kinematics for task-space	Tracking control requires adaptation of virtual trajectory
	Inverse dynamics required for torque-commanded robots**
Modular control architecture	Modular control architecture at the joint level	Modularity at the level of motion planning, interactive dynamics, and robot command
Inverse kinematics	Additional methods required for position- and torque-commanded robots	No inverse kinematics required
		Stable near (and even at) kinematic singularities
Managing kinematic redundancy	Must be explicitly handled	Redundancy managed by simple superposition of Mechanical Impedances
Robustness	Additional feedback control required for stable physical interaction	Inherently robust against passive environments
		Physical interaction explicitly regulated via Mechanical Impedances

*A map from trajectory generated by DMP to torque command is required (Section 3.1). **Methods for managing kinematic singularity (Section 3.3) and redundancy (Section 3.9, 3.10) must be included.

• • For DMPs, a low-gain PD controller is superimposed to manage uncertainty and physical contact (Section 3.1). EDAs include mechanical impedance as a separate primitive to manage physical interaction (Section 2.2.3). The dynamics of physical interaction can be controlled by modulating mechanical impedance (Section 3.4).

Once a detailed comparison has been made, we next show how DMPs and EDAs may be combined to leverage the best of both approaches (Section 4). Two major implementation examples on a KUKA LBR iiwa are provided.

2. Theory

In this Section, we provide an overview of DMPs and EDAs. For simplicity, we consider a system with a single DOF. A generalization to systems with multiple DOFs is presented in Section 3.

2.1. Dynamic movement primitives

DMPs, introduced by Schaal (1999, 2006); Ijspeert et al. (2013), parameterize the controller using the following three elements (Saveriano et al., 2023): a canonical system (Section 2.1.1), a nonlinear forcing term (Section 2.1.2), and a transformation system (Section 2.1.3). To generate discrete and rhythmic movements, two distinct definitions exist for the canonical system and nonlinear forcing term. For clarification, labels “Discrete” and “Rhythmic” are added next to the equation.

2.1.1. Canonical system

A canonical system $s:{\mathbb{R}}_{\ge 0}\to {\mathbb{R}}_{\ge 0}$ is a scalar variable governed by a first-order differential equation: $\tau \dot{s}\left(t\right)=\{\begin{array}{ll}-{\alpha }_{s}s\left(t\right)& \text{Discrete}\\ 1& \text{Rhythmic}\end{array}$ (1)In these equations, α _s is a positive constant, $t\in {\mathbb{R}}_{\ge 0}$ is time, and τ > 0 is a time constant.

For discrete movements, the canonical system is exponentially convergent to 0 with a closed-form solution s(t) = exp(−α _s t/τ)s(0). For rhythmic movements, the canonical system is a linear function of time s(t) = t/τ, but the modulo-2π operation is applied to ensure s(t) ∈ [0, 2π).

2.1.2. Nonlinear forcing term

A nonlinear forcing term $f:{\mathbb{R}}_{\ge 0}\to \mathbb{R}$ , which takes the canonical system s(t) as the function argument, is defined by: $f\left(s\left(t\right)\right)=\{\begin{array}{ll}\frac{{\sum }_{i=1}^N{w}_i{\varphi }_i\left(s\left(t\right)\right)}{{\sum }_{i=1}^N{\varphi }_i\left(s\left(t\right)\right)}s\left(t\right)\left(g-y_0\right)& \text{Discrete}\\ \frac{{\sum }_{i=1}^N{w}_i{\varphi }_i\left(s\left(t\right)\right)}{{\sum }_{i=1}^N{\varphi }_i\left(s\left(t\right)\right)}r& \text{Rhythmic}\end{array}$ (2)In equation (2), N is the number of basis functions; w _i is the weight of the ith basis function; y₀ and g are the initial and final positions of the discrete movement, respectively; r is the amplitude of the nonlinear forcing term for rhythmic movements; and ${\varphi }_i:{\mathbb{R}}_{\ge 0}\to \mathbb{R}$ is the ith basis function of the nonlinear forcing term: ${\varphi }_i\left(s\left(t\right)\right)=\{\begin{array}{ll}\exp \left\{-h_i{\left(s\left(t\right)-c_i\right)}^2\right\}& \text{Discrete}\\ \exp \left\{h_i\left(\cos \left(s\left(t\right)-c_i\right)-1\right)\right\}& \text{Rhythmic}\end{array}$ (3)In equation (3), the basis functions for discrete and rhythmic movements are Gaussian functions and von Mises functions, respectively (Ijspeert et al., 2013; Saveriano et al., 2023); c _i is the center of the ith basis function; and h _i is a positive constant that determines the width of the ith basis function.

The weights of the nonlinear forcing term w _i can be learned through various methods (Kober et al., 2013; Saveriano et al., 2023). For instance, Reinforcement Learning methods such as Natural Actor-Critic (Peters and Schaal, 2008) or PoWER (Muelling et al., 2010) can be used. Another prominent method is “Imitation Learning” (Ijspeert et al., 2013) (or “Learning from Demonstration”) (Billard et al., 2008; Calinon et al., 2007), where the weights are learned to follow (or imitate) a demonstrated trajectory. In this paper, we use Imitation Learning to learn the weights of the nonlinear forcing term. Further details on Imitation Learning are presented in Section 2.1.4.

2.1.3. Transformation system

The nonlinear forcing term f, with canonical system s as its function variable, is used as an input to the transformation system to generate trajectories with arbitrary complexity. In detail, a transformation system is a second-order differential equation, which is equivalent to a linear mass-spring-damper model with a nonlinear force input f(s(t)): $\begin{array}{l}\tau \dot{y}\left(t\right)=z\left(t\right)\\ \tau \dot{z}\left(t\right)={\alpha }_z\left\{{\beta }_z\left(g-y\left(t\right)\right)-z\left(t\right)\right\}+f\left(s\left(t\right)\right)\end{array}$ (4)In these equations, α _z and β _z are positive constants; y(t) and z(t) are state variables which correspond to position and (time-scaled) velocity of the transformation system, respectively; τ is the time constant used for the canonical system (equation (1)); for discrete movement, g is identical to equation (2); and for rhythmic movement, g is chosen to be the average of the rhythmic, repetitive movement (Ijspeert et al., 2013).

While any positive values of α _z and β _z can be used, usually the values of α _z and β _z are chosen such that if f(s(t)) is zero, the transformation system is critically damped (i.e., has repeated eigenvalues) for τ = 1 (i.e., β _z = α _z /4) (Ijspeert et al., 2013).

Using the canonical system s avoids the explicit time-dependency of the transformation system and results in an autonomous system (Ijspeert et al., 2013). Moreover, while both discrete and rhythmic movements are generated with the same transformation system, different choices of canonical systems s and nonlinear forcing terms f are made to produce those movements. Hence, to produce a combination of discrete and rhythmic movements, both discrete and rhythmic DMPs must be constructed (Section 3.7).

2.1.4. Imitation Learning

If the nonlinear forcing term is zero, that is, f(s(t)) = 0, the transformation system follows the response of a stable second-order linear system (equation (4)) with y(t) exponentially converging to g. Hence, y(t) can only represent a goal-directed discrete movement which follows the response of a stable linear system. To generate a wider range of trajectories (e.g., rhythmic movements), the weights of the nonlinear forcing terms must be learned.

For this, one prominent application of DMPs is “Imitation Learning,” also called “Learning from Demonstration” (Ijspeert et al., 2001, 2002, 2013; Schaal, 1999). Let y _des (t) be the desired trajectory that the robot aims to learn (or imitate). Then, y(t) = y _des (t) can be achieved by using the following nonlinear forcing term (equation (4)): $f_{target}\left(s\left(t\right)\right)={\tau }^2{\ddot{y}}_{des}\left(t\right)+{\alpha }_z\tau {\dot{y}}_{des}\left(t\right)+{\alpha }_z{\beta }_z\left(y_{des}\left(t\right)-g\right)$

If the analytic solution of y _des (t) and its derivatives are not known, Imitation Learning generates the whole continuous trajectory of y _des (t) from P sample points, (y _des (t _i ), ${\dot{y}}_{des}\left(t_i\right)$ , ${\ddot{y}}_{des}\left(t_i\right))$ for i ∈ {1, 2, …, P}, and the P sample points are used to find the best-fit weights of f(s(t)) (equation (3)) that matches f _target (s(t)).

The best-fit weight $w_i^{*}$ of the ith basis function is calculated using Locally Weighted regression (Atkeson et al., 1997; Ijspeert et al., 2013; Saveriano et al., 2023; Schaal and Atkeson, 1998): $w_i^{*}=\frac{{\mathbf{a}}^{\mathrm{T}}{\mathrm{\Phi }}_i{\mathbf{f}}_{target}}{{\mathbf{a}}^{\mathrm{T}}{\mathrm{\Phi }}_i\mathbf{a}}$ (5)where $\begin{array}{ll}\mathbf{a}& =\left[\begin{array}{l}{a}_1\\ a_2\\ ⋮\\ a_P\end{array}\right]\hspace{1em}{\mathbf{f}}_{target}=\left[\begin{array}{l}{f}_{target,1}\\ f_{target,2}\\ ⋮\\ f_{target,P}\end{array}\right]\\ {\mathrm{\Phi }}_i& =\left[\begin{array}{cccc}{\varphi }_i\left(s\left(t_1\right)\right)& & & 0\\ & {\varphi }_i\left(s\left(t_2\right)\right)& & \\ & & \ddots & \\ 0& & & {\varphi }_i\left(s\left(t_P\right)\right)\end{array}\right]\end{array}$

The elements of a and f _target are: $\begin{array}{ll}{a}_i\equiv a\left(t_i\right)& =\{\begin{array}{ll}s\left(t_i\right)\left(g-y_0\right)& \text{Discrete}\\ r& \text{Rhythmic}\end{array}\end{array}$ (6) $\begin{array}{ll}{f}_{target,i}& ={\tau }^2{\ddot{y}}_{des}\left(t_i\right)+{\alpha }_z\tau {\dot{y}}_{des}\left(t_i\right)+{\alpha }_z{\beta }_z\left(y_{des}\left(t_i\right)-g\right)\end{array}$

Along with Locally Weighted regression, one can also use Linear Least-Squares regression to find the best fit weights (Saveriano et al., 2019, 2023; Ude et al., 2014).

For Imitation Learning of discrete movements, the goal g and the initial position y₀ are set as g = y _des (t _P ) and y₀ = y _des (t₁), respectively; τ is chosen to be the duration of the discrete movement. For Imitation Learning of rhythmic movements, goal g is the midpoint of the minimum and maximum values of y _des (t₁), y _des (t₂), …, y _des (t _P ); τ is chosen to be the period of the demonstrated movement divided by 2π, hence the period of the rhythmic movement must be derived first (Ijspeert et al., 2013).

With these best-fit weights $w_i^{*}$ learned through Imitation Learning, the best-fit nonlinear forcing term f*(s(t)) is derived and used as the input to the transformation system to generate y _des (t), ${\dot{y}}_{des}\left(t\right)$ , and ${\ddot{y}}_{des}\left(t\right)$ .

One might wonder why spline methods are not used to derive y _des (t), ${\dot{y}}_{des}\left(t\right)$ , and ${\ddot{y}}_{des}\left(t\right)$ from the P sample points (Wada and Kawato, 2004). Splines are effective for smooth trajectory generation and are widely used in industrial robotics. However, spline methods do not allow online trajectory modulation of DMPs (Ijspeert et al., 2013) which is crucial to achieve multiple control tasks, for example, obstacle avoidance (Section 3.5) and sequencing discrete movements (Section 3.8). Compared to spline methods, Imitation Learning provides modularity for robot control. Once the best-fit weights are learned, these weights can be saved as a learned module, which can be reused to generate the learned trajectory (Section 3.3, 3.5). Finally, DMPs provide favorable invariance properties. Once the best-fit weights of the demonstrated trajectory are learned, the learned movement can be scaled in time (i.e., temporal invariance) and space (i.e., spatial invariance) without changing the qualitative property of the movement (Ijspeert et al., 2013).

If the P sample points of a demonstrated trajectory are provided, the best-fit weights which produce that demonstrated trajectory can be calculated with matrix algebra, referred to as “one-shot learning” (equation (5)). This process is called “batch regression” (Ijspeert et al., 2013), since all P data points should be collected to calculate the N weights.

The batch regression method assumes a predefined number of basis functions N and parameters c _i and h _i . While these parameters can be defined manually, the Locally Weighted Projection regression method (Vijayakumar and Schaal, 2000) can identify the necessary number of basis functions N, the center locations c _i , and width parameters h _i using nonparametric regression techniques.

Note that Imitation Learning for joint trajectories is easily scalable to high DOF systems (Atkeson et al., 2000; Ijspeert et al., 2002). For an n-DOF system, one can construct n transformation systems, each representing a joint trajectory. The n transformation systems are synchronized with a single canonical system (Ijspeert et al., 2013). With the Imitation Learning method, learning the best-fit weights is computationally efficient as it involves simple matrix algebra.

2.2. Elementary dynamic actions

EDAs, introduced by Hogan and Sternad (2012, 2013), consist of (at least) three distinct classes of primitives: submovements (Section 2.2.1) and oscillations (Section 2.2.2) as kinematic primitives, and mechanical impedances as interaction primitives (Section 2.2.3) (Figure 1).

Figure 1.

The three primitives of Elementary Dynamic Actions (EDAs). Submovements and oscillations correspond to kinematic primitives and mechanical impedances manage physical interaction.

2.2.1. Submovements

A submovement $x_0:{\mathbb{R}}_{\ge 0}\to \mathbb{R}$ is a smooth trajectory in which its time derivative is a unimodal function, that is, has a single peak value: ${\dot{x}}_0\left(t\right)=v\widehat{\sigma }\left(t\right)$ (7)In this equation, $\widehat{\sigma }:{\mathbb{R}}_{\ge 0}\to \mathbb{R}$ denotes a smooth unimodal basis function in which the integration over the whole time domain is 1, that is, ${\int }_0^{\infty }\widehat{\sigma }\left(t\right)dt=1$ ; and $v\mathbb{\in }\mathbb{R}$ is the velocity amplitude of the submovement.

Submovements model discrete motions, and therefore $\widehat{\sigma }\left(t\right)$ has a finite support, that is, there exists T > 0 such that $\widehat{\sigma }\left(t\right)=0$ for t ≥ T (Hogan and Sternad, 2007). The shape of $\widehat{\sigma }\left(t\right)$ can either be symmetric or not.

Given an initial condition, x₀(t = 0) ≡ x_0,i, the position of the submovement is defined by: $x_0\left(t\right)=x_{0,i}+v{f}_{\sigma }\left(t\right)$ (8)In this equation, $f_{\sigma }:{\mathbb{R}}_{\ge 0}\to {\mathbb{R}}_{\ge 0}$ is an integral of $\widehat{\sigma }\left(t\right)$ , that is, $f_{\sigma }\left(t\right)={\int }_0^t\widehat{\sigma }\left(s\right)ds$ ; accounting for the definition of $\widehat{\sigma }\left(t\right)$ , ∀t ≥ T : f _σ (t) = 1.

Note that x₀(t) serves as a general notation, encompassing various representations such as a joint-space trajectory (Section 3.2), or a trajectory for end-effector position (Section 3.3) and end-effector orientation (Section 3.10).

As discussed in Hogan and Sternad (2012), the definition of submovement is not only to provide a detailed mathematical formulation for practical application, but also as an account of observable motor behavior of humans. Prior studies have shown that planar reaching movements of unimpaired subjects followed a highly stereotyped unimodal speed profile (Atkeson and Hollerbach, 1985; Flash and Henis, 1991; Hogan and Flash, 1987; Park et al., 2017; Rohrer et al., 2002, 2004). The mathematical definition of submovements (equation (7)) accounts for these observable motor behaviors, while also enabling their use to generate goal-directed discrete movement.

2.2.2. Oscillations

An oscillation $x_0:{\mathbb{R}}_{\ge 0}\to \mathbb{R}$ is a smooth non-zero trajectory which is a periodic function: $\forall t>0:\exists T>0:x_0\left(t\right)=x_0\left(t+T\right)$

Note that this definition of oscillation can be too strict and the definition of oscillation can be expanded to almost-periodic functions (Hogan and Sternad, 2012). For our purposes, it is sufficient to think of an oscillation as a periodic function. Compared to submovements, oscillations model rhythmic and repetitive motions.

2.2.3. Mechanical impedances

Mechanical impedance Z is an operator which maps (generalized) displacement $\mathrm{\Delta }x\left(t\right)\mathbb{\in }\mathbb{R}$ to (generalized) force $F\left(t\right)\mathbb{\in }\mathbb{R}$ (Hogan, 2017; Hogan and Buerger, 2018): $Z:\mathrm{\Delta }x\left(t\right)\to F\left(t\right)$ (9)In detail, Δx(t) is the displacement of an actual trajectory of (generalized) position x(t) from a virtual trajectory x₀(t) to which the mechanical impedance is connected, that is, Δx(t) = x₀(t) − x(t). Loosely speaking, mechanical impedance is a generalization of stiffness to encompass nonlinear dynamic behavior (Hogan, 2017).

The impedance operator of equation (9) can denote both a map from joint displacement to torque or a map from end-effector (generalized) displacement to (generalized) force. The former impedance operator is often referred to as “joint-space impedance,” and the latter is often referred to as “task-space impedance.” For task-space impedance, both translational (Hogan, 1985) and rotational displacement (Caccavale et al., 1998) of the end-effector can be considered separately. The former is referred to as “position task-space impedance,” and the latter is referred to as “orientation task-space impedance.”

Along with the kinematic primitives, that is, submovements and oscillations, EDAs include mechanical impedance as a distinct primitive to manage physical interaction (Dietrich and Hogan, 2022; Hogan, 2017, 2022; Hogan and Buerger, 2018). The dynamics of physical interaction can be controlled by modulating mechanical impedance. For instance, tactile exploration and manipulation of fragile objects should evoke the use of low stiffness, while tasks such as drilling a hole on a surface require high stiffness for object stabilization (Hogan and Buerger, 2018).

Under the assumption that the environment is an admittance, mechanical impedances can be linearly superimposed even though each mechanical impedance is a nonlinear operator. This is the superposition principle of mechanical impedances (Hogan, 1985, 2017): $Z=\sum Z_i$ (10)

This principle provides a modular framework for robot control that can simplify multiple control tasks, for example, obstacle avoidance (Section 3.5) or managing kinematic redundancy (Section 3.9, 3.10).

Note that the impedance operators of equation (10) can include transformation maps. For instance, to superimpose a joint-space impedance and a task-space impedance at the torque level, the task-space impedance is multiplied by a Jacobian transpose to map from end-effector (generalized) force to joint torques (Section 3.1.2).

The choice of mechanical impedance decides whether the virtual trajectory x₀(t) is an attractor or repeller. The former can be used to produce discrete point-to-point movements (Section 3.2, 3.3), whereas the latter can be exploited for obstacle avoidance (Section 3.5) (Andrews and Hogan, 1983; Hjorth et al., 2020; Khatib, 1986; Newman, 1987).

2.2.4. Norton equivalent network model

The three distinct classes of EDAs—submovements, oscillations, and mechanical impedances—may be combined using a Norton equivalent network model (Hogan, 2017), which provides an effective framework to relate these classes of primitives (Figure 2).

Figure 2.

The three Elementary Dynamic Actions (EDAs) combined using a Norton equivalent network model. The virtual trajectory x₀(t) (yellow box) consists of submovements (orange box) and/or oscillations (blue box), and mechanical impedance Z (green box) regulates the interactive dynamics.

In detail, the forward-path dynamics (Figure 2) specifies virtual trajectory x₀(t), which consists of submovements and/or oscillations. The interactive dynamic, which consists of mechanical impedances Z, determines the generalized force output F(t) with the generalized displacement input Δx(t). Since the force output F(t) is determined by the choice of virtual trajectory x₀(t) and mechanical impedance Z, a key objective of EDAs is to find appropriate choices of x₀(t) and Z to produce the desired robot behavior.

Note that submovements and/or oscillations can be directly combined at the level of the virtual trajectory x₀(t), which thereby provides modularity at the kinematic level. As with the superposition principle of mechanical impedances (Section 2.2.3), the modular property at the kinematic level provides several benefits for multiple control tasks, for example, combining discrete and rhythmic movements (Section 3.7) and sequencing discrete movements (Section 3.8).

As shown in Figure 2, EDAs neither control x(t) (i.e., position control) nor F(t) (i.e., force/torque control) directly. Hence, EDAs are fundamentally different from position control (Arimoto, 1984; Tanner, 1981), force/torque control (Whitney, 1977), or hybrid position/force control methods (Raibert and Craig, 1981). This is one of the key ideas of impedance control. This is also the reason why we have chosen the terminology “virtual trajectory” for x₀(t) instead of a “desired” or a “reference trajectory.” Compared to tracking control methods which aim to follow a reference trajectory, x₀(t) of EDAs is simply a virtual trajectory to which the impedances are connected.

This property of EDAs has several benefits for robot control with physical interaction. Compared to x(t) and F(t) that depend on the environment or object with which the robot interacts, the virtual trajectory x₀(t) and impedance operator Z can be modulated “independently,” that is, regardless of the environment or the manipulated object (Hogan and Buerger, 2018). For instance, force/torque control cannot be used for free-space motions and position control cannot be used in contact with a kinematically constrained environment. EDAs can be used for both cases since neither x(t) nor F(t) is controlled directly.

Note that the Norton equivalent network model separates forward-path dynamics (virtual trajectory x₀(t)) from interactive dynamics (mechanical impedance Z). Hence, parallel optimization of x₀(t) and Z can be conducted. This has computational advantages for real-time control of robots with many DOF (Lachner et al., 2022).

3. Comparison of the two approaches

In this Section, a detailed comparison between DMPs and EDAs is presented. To emphasize the similarities and differences between the two approaches, multiple simulation examples using the MuJoCo physics engine (Version 1.50) (Todorov et al., 2012) are presented. The code is available at https://github.com/mosesnah-shared/DMP-comparison. ³ A list of examples, ordered in progressive complexity, is shown below:

• A goal-directed discrete movement in joint-space (Section 3.2).

• A goal-directed discrete movement for task-space position (Section 3.3).

• A goal-directed discrete movement for task-space position, with unexpected physical contact (Section 3.4).

• A goal-directed discrete movement for task-space position, including obstacle avoidance (Section 3.5).

• Rhythmic movement, both in joint-space and task-space position (Section 3.6).

• Combination of discrete and rhythmic movements, both in joint-space and task-space position (Section 3.7).

• A sequence of discrete movements for task-space position (Section 3.8).

• A single (or sequence of) discrete movement(s) for task-space position, while managing kinematic redundancy (Section 3.9).

• Discrete movement for task-space position and orientation, while managing kinematic redundancy (Section 3.10).

Some of the examples reproduce human-subject experiments in motor control research, for example, Burdet et al. (2001) for Section 3.3, and Flash and Henis (1991) for Section 3.8.

While, in general, position (or motion) control can be used to encode motor primitives, we will focus on the control of torque-actuated robots. Position control would create challenges that restrict the set of tasks that can be achieved (Hogan, 2022). For instance, one of the challenges is the kinematic transformation from task-space coordinates to the robot’s joint configuration, which complicates control in task-space. Another challenge is for tasks involving contact and physical interaction, which requires some level of compliance of the robotic manipulator. For the nine simulation examples, we highlight the challenges when using position-actuated robots, for example, managing contact and physical interaction (Section 3.4) and managing kinematic redundancy (Section 3.9, 3.10). We show that torque-actuated robots can address these tasks without imposing such challenges. Further discussion is deferred to Section 5.2.1.

Given a torque-actuated open-chain n-DOF robot manipulator, its dynamics is governed by the following differential equation (Murray et al., 1994): $\mathbf{M}\left(\mathbf{q}\right)\ddot{\mathbf{q}}+\mathbf{C}\left(\mathbf{q},\stackrel{.}{\mathbf{q}}\right)\stackrel{.}{\mathbf{q}}+\mathbf{G}\left(\mathbf{q}\right)={\mathit{\tau }}_{in}\left(t\right)+{\mathit{\tau }}_{ext}\left(t\right)$ (11)In this equation, $\mathbf{q}\equiv \mathbf{q}\left(t\right)\in {\mathbb{R}}^n$ is the joint trajectory of the robot; $\mathbf{M}\left(\mathbf{q}\right)\in {\mathbb{R}}^{n\times n}$ and $\mathbf{C}\left(\mathbf{q},\stackrel{.}{\mathbf{q}}\right)\in {\mathbb{R}}^{n\times n}$ are the mass and centrifugal/Coriolis matrices, respectively; $\mathbf{G}\left(\mathbf{q}\right)\in {\mathbb{R}}^n$ is the vector arising through gravitational potential energy; ${\mathit{\tau }}_{in}\left(t\right)\in {\mathbb{R}}^n$ is the torque input; ${\mathit{\tau }}_{ext}\left(t\right)\in {\mathbb{R}}^n$ is the resultant effect of external forces expressed as torque. For torque-controlled robots, the goal is to determine the torque input τ _in (t) which produces desired robot behavior. For brevity and to avoid clutter, we often omit argument t.

In this paper, we consider controlling both position and orientation of the end-effector. For end-effector position, we use $\mathbf{p}\equiv \mathbf{p}\left(t\right)\in {\mathbb{R}}^3$ to denote the three-dimensional Cartesian position of the end-effector and h _p to represent its Forward Kinematic Map, that is, p = h _p (q). For end-effector orientation, we use R ≡ R(t) ∈ SO(3) to denote the spatial orientation of the end-effector, where SO(3) denotes the Special Orthogonal Group in three-dimensional space (Lynch and Park, 2017; Murray et al., 1994). To represent the Forward Kinematic Map of the robot for orientation, we use h _r , that is, R = h _r (q).

Note that one can also represent spatial orientation with unit quaternions $\overrightarrow{\mathbf{q}}\in {\mathbb{H}}_1$ , ⁴ where ${\mathbb{H}}_1$ denotes the set of quaternions with unit norm (Caccavale et al., 1998; Shuster, 1993). In this case, we denote the end-effector orientation with a quadruple $\overrightarrow{\mathbf{q}}=\left(\eta ,\mathit{ϵ}\right)\in {\mathbb{H}}_1$ , where $\eta \mathbb{\in }\mathbb{R}$ (respectively $\mathit{ϵ}\in {\mathbb{R}}^3$ ) is the real (respectively imaginary) part of the unit quaternion (Caccavale et al., 1998; Lachner et al., 2022). One can convert between R and unit quaternion, using numerically stable methods presented in Shepperd (1978) and Allmendinger (2015).

Except for tasks with physical contact (e.g., Section 3.4), we assume τ _ext (t) = 0. Finally, we assume that gravitational force G(q) is compensated by the controller and can be neglected (equation (11)).

3.1. The existence of an inverse dynamics model

We first show that for torque-actuated robots, DMPs require an inverse dynamics model, whereas EDAs do not.

3.1.1. Dynamic movement primitives

For DMPs, the transformation system (equation (4)) represents kinematic relations. Hence for a torque-actuated robot, the approach requires an inverse dynamics model, which determines the feedforward joint torques τ _in (t) to generate the desired joint trajectory specified by $\ddot{\mathbf{q}}\left(t\right)$ , $\stackrel{.}{\mathbf{q}}\left(t\right)$ , and q(t) (Ijspeert et al., 2013; Saveriano et al., 2023). The inverse dynamics model requires an exact model of the robot manipulator, that is, exact values of M(q) and $\mathbf{C}\left(\mathbf{q},\stackrel{.}{\mathbf{q}}\right)$ .

The requirement of an inverse dynamics model also implies that the DMP approach is in principle, a nonreactive feedforward control approach. To be robust against uncertainty or unexpected physical contact, a low-gain feedback control (e.g., PD control) is added to the feedforward inverse dynamics controller (Pastor et al., 2013; Schaal et al., 2007) (Section 3.3, 3.4). Moreover, for control in task-space where the maps from p(t), R(t), and their higher derivatives to $\ddot{\mathbf{q}}\left(t\right)$ , $\stackrel{.}{\mathbf{q}}\left(t\right)$ , and q(t) are not always well defined, additional control methods must be employed. For instance, to solve inverse kinematics, methods presented by Whitney (1969); Liegeois (1977); Maciejewski and Klein (1988) are used (Nakanishi et al., 2008; Park et al., 2008; Pastor et al., 2009) (Section 3.9.3, 3.10.3). For tracking control in task-space, feedback control methods such as sliding-mode control (Nakanishi et al., 2008; Slotine and Li, 1991) or its modification (Ceccarelli, 2008) are employed (Section 3.9.3, 3.10.3).

3.1.2. Elementary dynamic actions

For EDAs, an inverse dynamics model is not required for a torque-actuated robot. Hence, an exact model of the robot, that is, exact computation of M(q) and $\mathbf{C}\left(\mathbf{q},\stackrel{.}{\mathbf{q}}\right)$ is not essential, unless advanced applications such as inertia shaping (Dietrich and Ott, 2019; Khatib, 1995) or regulating the physical interaction with the environment (Lachner et al., 2021) (Section 3.4.2) are desired. Instead, the input torque command τ _in (t) is determined by superimposing mechanical impedances (Section 2.2.3): $\begin{array}{ll}{\mathit{\tau }}_{in}\left(t\right)& ={\displaystyle \sum }_i{\mathbf{Z}}_{q,i}\left(\mathrm{\Delta }\mathbf{q}\left(t\right)\right)+{\displaystyle \sum }_i{\mathbf{J}}_p^{\mathrm{T}}\left(\mathbf{q}\left(t\right)\right){\mathbf{Z}}_{p,i}\left(\mathrm{\Delta }\mathbf{p}\left(t\right),t\right)\\ & +{\displaystyle \sum }_i{\mathbf{J}}_r^{\mathrm{T}}\left(\mathbf{q}\left(t\right)\right){\mathbf{Z}}_{r,i}\left(\mathrm{\Delta }\mathbf{R}\left(t\right),t\right)\end{array}$ In this equation, ${\mathbf{J}}_p\left(\mathbf{q}\left(t\right)\right)\equiv {\mathbf{J}}_p\left(\mathbf{q}\right)\in {\mathbb{R}}^{3\times n}$ and ${\mathbf{J}}_r\left(\mathbf{q}\left(t\right)\right)\equiv {\mathbf{J}}_r\left(\mathbf{q}\right)\in {\mathbb{R}}^{3\times n}$ are the Jacobian matrices for task-space position and orientation, respectively (Siciliano et al., 2008; Siciliano and Khatib, 2016); the former is often referred to as the Analytical Jacobian matrix, which can be derived by the partial derivatives of h _p (q) with respect to q, that is, J _p (q) ≡ $\frac{\partial {\mathbf{h}}_p}{\partial \mathbf{q}}\left(\mathbf{q}\right)$ ; the latter is referred to as the Geometric Jacobian matrix (Lynch and Park, 2017; Murray et al., 1994; Siciliano et al., 2008; Siciliano and Khatib, 2016); Z _q , Z _p , and Z _r denote joint-space, position task-space, and orientation task-space mechanical impedances, respectively.

Note that Δq(t) and Δp(t) denote a simple subtraction between virtual and actual trajectories, that is, Δq(t) ≡ q₀(t) − q(t) and Δp(t) ≡ p₀(t) − p(t), where q₀(t) and p₀(t) denote the virtual trajectories for joint-space and task-space position, respectively. For task-space orientation, ΔR(t) ≡ R^T(t)R₀(t), where R₀(t) denotes the virtual trajectory for task-space orientation.

Compared to DMPs, the EDA approach is in principle, a reactive feedback control method. Instead of an inverse dynamics model, measurements q(t) and a Forward Kinematics Map of the robot are required. Given these values, one can react to the environment by modulating impedance and/or the virtual trajectories to regulate the dynamics of physical interaction (Section 3.4). Moreover, with appropriate choices of Z _r , Z _p , Z _q and R₀(t), p₀(t), q₀(t), the controller preserves passivity, which thereby provides robustness against uncertainty and external disturbances (Section 3.4).

3.2. A goal-directed discrete movement in joint-space

We consider designing a controller to generate a goal-directed discrete movement planned in joint-space coordinates.

3.2.1. Dynamic movement primitives

The movement of each joint is represented by a transformation system. The n transformation systems are synchronized by using a single discrete canonical system as the input to the n nonlinear forcing terms (Ijspeert et al., 2013) (Section 2.1.4): $\begin{array}{ll}{\tau }^2\ddot{\mathbf{q}}\left(t\right)& =-{\alpha }_z{\beta }_z\left\{\mathbf{q}\left(t\right)-\mathbf{g}\right\}-{\alpha }_z\tau \stackrel{.}{\mathbf{q}}\left(t\right)+\mathbf{f}\left(s\left(t\right)\right)\\ \tau \dot{s}\left(t\right)& =-{\alpha }_{s}s\left(t\right)\end{array}$ (12)

Note that different values of α _z and β _z can be used for each transformation system. Nevertheless, it is sufficient to use identical values of α _z and β _z for n transformation systems.

Without using a nonlinear forcing term, that is, f(s(t)) = 0, the goal-directed discrete movement can be generated by setting g as the goal location. As a result, q(t) generates a motion of a stable second-order linear system which converges to g.

In case we want to generate a goal-directed discrete movement that also follows a specific joint trajectory, Imitation Learning can be used. Let q _des (t) be the desired discrete joint-trajectory with duration T, such that q _des (T) = g and τ = T. The best-fit force f*(s(t)) is calculated and used as the input to n transformation systems to produce ${\ddot{\mathbf{q}}}_{des}\left(t\right)$ (equation (12)). With the initial conditions ${\mathbf{q}}_{des}\left(0\right),{\stackrel{.}{\mathbf{q}}}_{des}\left(0\right)$ , the trajectories of q _des (t), ${\stackrel{.}{\mathbf{q}}}_{des}\left(t\right)$ are calculated using numerical integration. The calculated q _des (t), ${\stackrel{.}{\mathbf{q}}}_{des}\left(t\right)$ ,and ${\ddot{\mathbf{q}}}_{des}\left(t\right)$ are the input to the inverse dynamics model to generate the necessary torque input τ _in (t).

3.2.2. Elementary dynamic actions

To generate a goal-directed discrete movement planned in joint-space coordinates, we construct the following controller: ${\mathit{\tau }}_{in}\left(t\right)={\mathbf{K}}_q\left\{{\mathbf{q}}_0\left(t\right)-\mathbf{q}\right\}+{\mathbf{B}}_q\left\{{\stackrel{.}{\mathbf{q}}}_0\left(t\right)-\stackrel{.}{\mathbf{q}}\right\}$ (13)In this equation, ${\mathbf{K}}_q,{\mathbf{B}}_q\in {\mathbb{R}}^{n\times n}$ are symmetric positive definite matrices which correspond to the joint stiffness and damping, respectively. This controller is a first-order joint-space impedance controller (Nah et al., 2020, 2021).

To generate a goal directed discrete movement, q₀(t) is chosen to be a submovement which ends at goal location g, that is, if the duration of submovement is T, then q₀(T) = g. For constant positive definite K _q , B _q matrices, and q(t) asymptotically converges to g (Slotine and Li, 1991; Takegaki, 1981).

3.2.3. Simulation example

Consider a 2-DOF planar robot model, where each link consists of a single uniform slender bar with mass and length of 1 kg and 1 m, respectively. Let the initial joint configuration be $\mathbf{q}\left(0\right)={\mathbf{q}}_i\in {\mathbb{R}}^2$ . The code script for this simulation is main_joint_discrete.py.

For q _des (t) of DMPs and q₀(t) of EDAs, both trajectories were set to be a minimum-jerk trajectory (Flash and Hogan, 1985; Hogan and Flash, 1987): $\begin{array}{l}{\mathbf{q}}_{des}\left(t\right),{\mathbf{q}}_0\left(t\right)=\{\begin{array}{ll}{\mathbf{q}}_i+\left(\mathbf{g}-{\mathbf{q}}_i\right)f_{MJT}\left(t\right)& 0\le t<T\\ \mathbf{g}& T\le t\end{array}\\ f_{MJT}\left(t\right)=10{\left(\frac{t}{T}\right)}^3-15{\left(\frac{t}{T}\right)}^4+6{\left(\frac{t}{T}\right)}^5\end{array}$ (14)

For submovement q₀(t), the parameters were x_0,i = q _i , v = g − q _i , and f _σ (t) = f _MJT (t) for 0 ≤ t ≤ T (equation (8)).

The simulation results are shown in Figure 3. DMPs generated the goal-directed discrete movement while achieving perfect tracking of the minimum-jerk trajectory. Once the weights of the nonlinear forcing terms are learned using Imitation Learning, one can regenerate the minimum-jerk trajectory by simply retrieving these learned weights. While the presented example used a minimum-jerk trajectory, Imitation Learning can be used to achieve tracking control of a trajectory with arbitrary complexity.

Figure 3.

A goal-directed discrete movement in joint-space for Dynamic Movement Primitives (DMPs, blue) and Elementary Dynamic Actions (EDAs, orange) (Section 3.2.3). The black dotted line is a minimum-jerk trajectory (equation (14)). Goal location: g = [1, 1] rad. Parameters of minimum-jerk trajectory: q _i = [0, 0] rad, T = 1s. Parameters of DMPs: α _z = 10, β _z = 2.5, α _s = 1.0, τ = T, N = 50, p = 100, c _i = exp(−α _s (i − 1)/(N − 1)) for i ∈ [1, 2, …, N], $h_i=1/{\left(c_{i+1}-c_i\right)}^2$ for i ∈ [1, 2, …, N − 1], h _N = h_N−1. Parameters of EDAs: K _q = 150I₂ N⋅m/rad, B _q = 50I₂ N⋅m⋅s/rad, where ${\mathbf{I}}_2\in {\mathbb{R}}^{2\times 2}$ is an identity matrix. For DMP, perfect tracking was achieved. For EDA, a non-negligible tracking error was observed.

In principle, perfect tracking may be achieved using DMPs. However, in practice, one must compensate for the inaccuracy of the inverse dynamics model by superimposing a joint-space Proportional Derivative (PD) feedback controller. Further details on this method are deferred to Section 3.4.

For EDAs, a non-zero error between q₀(t) and q(t) was observed. Hence, perfect tracking of the minimum-jerk trajectory was not achieved. Nevertheless, EDAs enabled asymptotic convergence of q(t) toward the goal configuration g without the need for an inverse dynamics model. Moreover, one can reduce the tracking error in joint-space by increasing the joint stiffness and damping values (Figure 4).

Figure 4.

Joint-space tracking performance of EDAs with respect to different values of mechanical impedances. (1st, 2nd, and 3rd Columns) Performance with respect to different joint stiffness matrices K _q (equation (13)). Joint damping value: B _q = 50I₂ N⋅m⋅s/rad. Higher stiffness minimizes tracking error, while lower stiffness increases settling time and rise time. (4th, 5th, 6th Columns) Performance with respect to different joint damping matrices B _q (equation (13)). Joint stiffness value: K _q = 150I₂ N⋅m/rad. Higher damping minimizes tracking error, while lower damping increases overshoot and emphasizes oscillation. For the 2nd and 5th columns, orange lines and dotted black lines are identical to those in Figure 3.

3.3. A goal-directed discrete movement for task-space position

We next design a controller to generate a goal-directed discrete movement of the end-effector’s position. For this Section, we assume no kinematic redundancy of the robot model, that is, the Jacobian matrix J _p (q) is a square matrix. Moreover, we only consider the end-effector’s position for control. A method to manage kinematic redundancy is considered in Section 3.9. The control of the end-effector’s orientation is considered in Section 3.10. One of the goal positions is located at a kinematic singularity, that is, a fully stretched configuration.

3.3.1. Dynamic movement primitives

For DMPs, the task can be achieved by representing the end-effector trajectory p(t) with the transformation system (Pastor et al., 2009): $\begin{array}{l}{\tau }^2\ddot{\mathbf{p}}\left(t\right)=-{\alpha }_z{\beta }_z\left(\mathbf{p}\left(t\right)-\mathbf{g}\right)-{\alpha }_z\tau \stackrel{.}{\mathbf{p}}\left(t\right)+\mathbf{f}\left(s\left(t\right)\right)\\ \tau \dot{s}\left(t\right)=-{\alpha }_{s}s\left(t\right)\end{array}$ (15)

Once the desired end-effector trajectories p _des , ${\stackrel{.}{\mathbf{p}}}_{des}$ , and ${\ddot{\mathbf{p}}}_{des}$ are computed with or without Imitation Learning (as discussed in Section 3.2), the inverse kinematics are used to calculate the corresponding joint trajectories: $\begin{array}{l}{\mathbf{q}}_{des}={\mathbf{h}}_p^{-1}\left({\mathbf{p}}_{des}\right)\\ {\stackrel{.}{\mathbf{q}}}_{des}={\mathbf{J}}_p^{-1}\left({\mathbf{q}}_{des}\right){\stackrel{.}{\mathbf{p}}}_{des}\\ {\ddot{\mathbf{q}}}_{des}={\mathbf{J}}_p^{-1}\left({\mathbf{q}}_{des}\right)\left\{{\ddot{\mathbf{p}}}_{des}-{\stackrel{.}{\mathbf{J}}}_p\left({\mathbf{q}}_{des}\right){\stackrel{.}{\mathbf{q}}}_{des}\right\}\end{array}$

The calculated ${\ddot{\mathbf{q}}}_{des},{\stackrel{.}{\mathbf{q}}}_{des},{\mathbf{q}}_{des}$ are used as the input to the inverse dynamics model to calculate τ _in (t) (Section 3.1).

Note that the presented method for inverse kinematics cannot be used for a kinematically redundant robot (with fewer task-space than joint-space DOFs). To manage kinematic redundancy, along with the feedforward torque command from the inverse dynamics model, an additional feedback controller should be employed (Nakanishi et al., 2008; Pastor et al., 2009; Slotine and Li, 1991). This is considered further in Section 3.9 and Section 3.10. Moreover, to handle kinematic singularity, that is, when J _p (q) loses rank and its inverse is not defined, methods such as Damped Least-squares inverse must be employed in addition (Buss and Kim, 2005; Chiaverini et al., 1994; Nakamura and Hanafusa, 1986).

3.3.2. Elementary dynamic actions

To generate a goal-directed discrete movement planned in task-space coordinates, we construct the following controller: ${\mathit{\tau }}_{in}\left(t\right)={\mathbf{J}}_p^{\mathrm{T}}\left(\mathbf{q}\right)\left[{\mathbf{K}}_p\left\{{\mathbf{p}}_0\left(t\right)-\mathbf{p}\right\}+{\mathbf{B}}_p\left\{{\stackrel{.}{\mathbf{p}}}_0\left(t\right)-\stackrel{.}{\mathbf{p}}\right\}\right]$ (16)In this equation, ${\mathbf{K}}_p,{\mathbf{B}}_p\in {\mathbb{R}}^{3\times 3}$ ( ${\mathbb{R}}^{2\times 2}$ for the planar simulations presented below) are constant symmetric positive definite matrices which correspond to translational stiffness and damping, respectively. This controller is a first-order impedance controller for task-space position (Hermus et al., 2021; Verdi, 2019).

Compared to DMPs (Section 3.3.1), this controller does not require a Jacobian inverse. Hence, the torque input is always well-defined near and even at kinematic singularities.

To generate a goal-directed discrete movement, as with the first-order joint-space impedance controller (Section 3.2), p₀(t) is chosen to be a submovement which ends at goal location g.

For constant positive definite K _p , B _p matrices, p(t) asymptotically converges to g (Takegaki, 1981). This also implies an asymptotic convergence of $\mathbf{q}\left(t\right)\to {\mathbf{h}}_p^{-1}\left(\mathbf{g}\right)$ .

3.3.3. Simulation example

The simulation example in Figure 5 reproduced the movement of the experiment conducted in Burdet et al. (2001). The goal-directed discrete movement was made in a direction away from the robot base, along the positive Y-axis direction (Figure 5). The amplitude of the movement was varied such that one of the movements reached a fully stretched configuration, that is, a configuration at kinematic singularity. The code script for this simulation is main_task_discrete.py.

Figure 5.

A goal-directed discrete movement for task-space position, for (A, C, E, F) Dynamic Movement Primitives (DMPs, blue) and (B, D, F) Elementary Dynamic Actions (EDAs, orange) (Section 3.3.3). (A, B, C, D, E) Goal-directed discrete movements in a direction along the positive Y-axis direction, but (C, D, E) reached kinematic singularity. The black dotted line is a minimum-jerk trajectory (equation (17)). (E) DMPs with Damped Least-squares inverse to manage kinematic singularity. As shown in the Figure, the upper-limb is fully stretched and completely overlaps the dotted black line. (F) Time t versus σ _min ((J _p (q)))/σ _max ((J _p (q))) for (C, D, E), where σ is a singular value of the matrix. For DMPs, the MuJoCo simulation halted near 0.99s. Parameters of a minimum-jerk trajectory: (A, B) p _i = [0.0, 0.52]m, g = [0.0, 1.72]m, T = 1.0s (C, D) p _i = [0.0, 0.52]m, g = [0.0, 2.00]m, T = 1.0s, where p _i is the initial end-effector position. Parameters of DMPs are identical to those in Figure 3. Parameters of EDAs: K _p = 60I₂ N/m, B _p = 20I₂ N⋅s/m. For DMP, perfect tracking was achieved when the motion was not near a kinematic singularity. However for (C, F) without Damped Least-squares inverse, when the goal g was at a kinematic singularity, the approach failed to achieve the task. (E, F) With Damped Least-squares inverse, the upper-limb remained stable near the kinematic singularity and the fully-stretched configuration was achieved. For EDA, a non-negligible tracking error was observed, but goal-directed discrete movements were achieved for both (B) and (D). The approach based on EDA was not only stable near and at a kinematic singularity but even passed through it (F, near t = 1.0 and t = 2.4, depicted as grey circle markers).

We used the 2-DOF planar robot model from Section 3.2 that was constrained to move within the XY plane. For DMPs and EDAs, both p _des (t) and p₀(t) were chosen to be a minimum-jerk trajectory (equation (14)): $\begin{array}{l}{\mathbf{p}}_{des}\left(t\right),{\mathbf{p}}_0\left(t\right)=\{\begin{array}{ll}{\mathbf{p}}_i+\left(\mathbf{g}-{\mathbf{p}}_i\right)f_{MJT}\left(t\right)& 0\le t<T\\ \mathbf{g}& T\le t\end{array}\\ f_{MJT}\left(t\right)=10{\left(\frac{t}{T}\right)}^3-15{\left(\frac{t}{T}\right)}^4+6{\left(\frac{t}{T}\right)}^5\end{array}$ (17)In this equation, p _i is the initial position of the end-effector; for submovement p₀(t), the parameters are x_0,i = p _i , v = g − p _i , and f _σ (t) = f _MJT (t) for 0 ≤ t ≤ T (equation (8)).

As shown in Figure 5(A) and (B), both approaches successfully generated discrete movements that converged to the desired goal location. As discussed in Section 3.2.3, DMPs achieved the goal-directed discrete movement with perfect tracking. However, as in Section 3.2.3, perfect tracking requires an accurate model of the robot’s kinematics and dynamics. Hence in practice, one must compensate for the inaccuracy of these models by superimposing a joint-space PD feedback controller (Section 3.4).

For EDAs, goal-directed discrete movement was achieved but a tracking error was observed. However, as in Section 3.2.3 (Figure 4), one can reduce the tracking error in task-space by increasing the translational stiffness and damping values (Figure 6).

Figure 6.

Task-space position tracking performance of EDAs with respect to different values of mechanical impedance. (1st, 2nd, 3rd Columns) Performance with respect to different translational stiffness matrices K _p (equation (16)). Translational damping value: B _p = 20I₂ N⋅s/m. Higher stiffness reduces tracking error, while lower stiffness increases settling time and rise time. (Right Column) Performance with respect to different translational damping matrices B _p (equation (16)). Translational stiffness value: K _p = 60I₂ N/m. Higher damping reduces tracking error, while lower damping increases overshoot and emphasizes oscillation. For the 2nd and 5th columns, orange lines and dotted black lines are identical to those in Figure 5(B).

Note that EDAs did not require the Jacobian inverse. The benefit of this property was emphasized when the planar robot model reached for a fully stretched configuration. As shown in Figures 5(C) and 5(F), DMPs without Damped Least-squares inverse became numerically unstable when the robot model approached a kinematic singularity. Hence, as shown in Figures 5(E) and (F), Damped Least-squares inverse must be employed to remain stable near kinematic singularity. For EDAs, not only was the approach stable, but the approach even “passed-through” the kinematic singularity (Figure 5(F), near t = 1.0 and t = 2.4, depicted as grey circle markers) without any numerical instability of the controller.

Note that in Figure 5(D), using EDAs, the robot model oscillated back and forth between its “left-hand” and “right-hand” configurations, passing through the singularity multiple times. This occurred because at the singular configuration the controller included no effective damping, hence no means to dissipate the angular momentum of the robot links. This oscillation may be suppressed by adding a non-zero joint-space damping, an example of impedance superposition (Section 2.2.3, 3.5, 3.9). Moreover, with EDAs, one can “exploit” rather “avoid” kinematic singularity. For instance, if a task can be facilitated by using the left-hand configuration, EDAs enable “shifting” from right- to left-hand configuration.

3.4. Tasks with unexpected physical contact

We next consider a case where unexpected physical contact is made while conducting a point-to-point discrete reaching movement presented in Section 3.3.3 (Figure 5(A) and (B)).

3.4.1. Dynamic movement primitives

As discussed in Section 3.1, DMPs use a feedforward torque command calculated from the inverse dynamics model. To handle model uncertainties and possible instabilities for control tasks involving unexpected contact, a feedback torque command is superimposed on the feedforward torque command.

Let τ _ff (t) be the feedforward torque command (Section 3.1, 3.3) to produce p _des (t). A feedback torque command τ _fb (t) based on a low-gain joint-space PD controller was superimposed on τ _ff (t): ${\mathit{\tau }}_{fb}\left(t\right)={\mathbf{K}}_q\left\{{\mathbf{q}}_{des}\left(t\right)-\mathbf{q}\right\}+{\mathbf{B}}_q\left\{{\stackrel{.}{\mathbf{q}}}_{des}\left(t\right)-\stackrel{.}{\mathbf{q}}\right\}$ (18) ${\mathit{\tau }}_{in}\left(t\right)={\mathit{\tau }}_{ff}\left(t\right)+{\mathit{\tau }}_{fb}\left(t\right)$ (19)

The gain values for K _q and B _q were manually chosen to be sufficiently small to manage unexpected contacts (Schaal et al., 2007). Moreover, ${\mathbf{q}}_{des}\left(t\right),{\stackrel{.}{\mathbf{q}}}_{des}\left(t\right)$ were derived by inverse kinematics of p _des (t), ${\stackrel{.}{\mathbf{p}}}_{des}\left(t\right)$ (Section 3.3.1). The gains K _q and B _q can also be determined by stochastic optimal control (Buchli et al., 2011; Theodorou et al., 2010).

Note that for ideal torque-actuated robots, the joint-space PD controller (equation (18)) is identical to a first-order joint-space impedance controller of EDAs (equation (13)). However, it would be a mistake to conclude that PD control is identical to impedance control (Won et al., 1997). Further discussion is deferred to Section 5.2.3.

3.4.2. Elementary dynamic actions

As discussed in Hogan (1985, 2022), an impedance controller is robust against unexpected physical contact with passive environments. While it is common to use constant mechanical impedances, mechanical impedance can be modulated to regulate the dynamics of physical interaction (Lachner et al., 2021).

In detail, the first-order position task-space impedance controller of equation (16) can be adapted by modulating the translational stiffness and damping values: ${\mathit{\tau }}_{in}={\mathbf{J}}_p^{\mathrm{T}}\left(\mathbf{q}\right)\left[{\mathbf{K}}_p^{\prime }\left(\lambda \right)\left\{{\mathbf{p}}_0-\mathbf{p}\right\}+{\mathbf{B}}_p^{\prime }\left(\lambda \right)\left\{{\stackrel{.}{\mathbf{p}}}_0-\stackrel{.}{\mathbf{p}}\right\}\right]$ (20)In this equation, ${\mathbf{K}}_p^{\prime }\left(\lambda \right)=\lambda {\mathbf{K}}_p$ and ${\mathbf{B}}_p^{\prime }\left(\lambda \right)=c\lambda {\mathbf{K}}_p$ ; K _p is a constant symmetric positive definite matrix; λ ∈ [0, 1], and c is a positive constant that determines the damping ratio.

With a slight abuse of notation, ${\mathbf{K}}_p^{\prime }\left(\lambda \right)$ and ${\mathbf{B}}_p^{\prime }\left(\lambda \right)$ were modulated via $\lambda \equiv \lambda \left(\mathcal{T},\mathcal{U},{\mathcal{L}}_{max}\right)$ , which is a function of the kinetic energy of the robot $\mathcal{T}\equiv \mathcal{T}\left(\mathbf{q}\left(t\right),\stackrel{.}{\mathbf{q}}\left(t\right)\right)$ , elastic potential energy $\mathcal{U}\equiv \mathcal{U}\left(\mathrm{\Delta }\mathbf{p}\left(t\right),{\mathbf{K}}_p\right)$ due to the translational stiffness K _p and displacement Δp(t), and the energy threshold ${\mathcal{L}}_{max}$ : $\lambda =\{\begin{array}{cc}1& \text{if}{\mathcal{L}}_c\left(t,\lambda \right)\le {\mathcal{L}}_{\max }\\ \max \left(\frac{1}{\mathcal{U}}\left({\mathcal{L}}_{\max }-\mathcal{T}\right),0\right)& \text{if}{\mathcal{L}}_c\left(t,\lambda \right)>{\mathcal{L}}_{\max }\end{array}$ where: $\mathcal{T}\left(\mathbf{q},\stackrel{.}{\mathbf{q}}\right)=\frac{1}{2}{\stackrel{.}{\mathbf{q}}}^{\mathrm{T}}\mathbf{M}\left(\mathbf{q}\right)\stackrel{.}{\mathbf{q}}\mathcal{U}\left(\mathrm{\Delta }\mathbf{p},{\mathbf{K}}_p\right)=\frac{1}{2}\mathrm{\Delta }{\mathbf{p}}^{\mathrm{T}}{\mathbf{K}}_p\mathrm{\Delta }\mathbf{p}{\mathcal{L}}_c\left(t,\lambda \right)=\mathcal{T}\left(\mathbf{q},\stackrel{.}{\mathbf{q}}\right)+\mathcal{U}\left(\mathrm{\Delta }\mathbf{p}\right)$

This controller limits the impact of an unexpected contact, for example, during physical Human-Robot Interaction (pHRI) (Lachner, 2022; Lachner et al., 2021), by bounding the total energy of the robot via λ. For $\mathcal{T}\le {\mathcal{L}}_{\max }$ , the total energy of the robot is always less than or equal to ${\mathcal{L}}_{\max }$ . To avoid negative stiffness and damping matrices emerging from the condition $\mathcal{T}>{\mathcal{L}}_{\max }$ , λ is set to be zero via the max-function.

For pHRI, the value ${\mathcal{L}}_{\max }$ is defined by standards and regulations (International Organization for Standardization, 2016) and is dependent on the human body part that is in contact with the robot.

While the presented controller is a simplified example, a more advanced application exists where the damping term can be modulated as a function of stiffness and inertia (Albu-Schaffer et al., 2003). Moreover, the dissipative behavior of the robot can be modulated to limit the robot power, for example, by using “damping injection” (Stramigioli, 2015).

Note that for this advanced application, EDAs require the mass matrix of the robot M(q) for the calculation of λ which regulates the translational stiffness ${\mathbf{K}}_p^{\prime }\left(\lambda \right)$ and damping matrices ${\mathbf{B}}_p^{\prime }\left(\lambda \right)$ (equation (20)).

3.4.3. Simulation example

As in Section 3.3.3, we used a 2-DOF planar robot model to generate a goal-directed discrete movement in task-space coordinates. In this example, a square-shaped obstacle was placed to block the robot path. A few seconds after the first contact with the robot, the obstacle was moved aside and the robot could continue its motion (Andrews and Hogan, 1983; Newman, 1987). The code script for this simulation is main_unexpected_contact.py.

For DMPs, without a low-gain PD feedback controller, the learned weights from Section 3.3.3 were reused. For this controller, the robot model bounced back from the obstacle due to contact and failed to reach the goal (Figure 7(B)). Hence, it was necessary to add a low-gain PD controller (Figure 7(C)). This example demonstrated the modular property of DMPs, since the learned feedforward torque controller from Section 3.3 was reused without modification, and superimposed upon an additional feedback controller.

Figure 7.

A goal-directed discrete movement with unexpected physical contact for (A, B, C, D) Dynamic Movement Primitives (DMPs, blue) and (E, F, G, H) Elementary Dynamic Actions (EDAs, orange) (Section 3.4.3). For DMP: (A) Moment of contact. (A → B) DMP without PD controller (A → C) DMP with PD controller. (D) Time versus Y-coordinate of the end-effector for A → B (dotted line) and A → C (filled line). For EDA: (E) Moment of contact. (E → F) EDA without impedance modulation (E → G) EDA with impedance modulation. (H) Time versus Y-coordinate of the end-effector for E → F (dotted line) and E → G (filled line). With impedance modulation of E → G, a slower movement than E → F was achieved. (D, H) The black dotted lines are a minimum-jerk trajectory (equation (17)). Parameters of the minimum-jerk trajectory, DMPs and EDAs are identical to those in Figure 5. Parameters of the PD controller of DMPs (equation (18)): K _q = 50I₂ N⋅m/rad, B _q = 30I₂ N⋅m⋅s/rad. Smaller gain values than Figure 3 were chosen to achieve a relatively low-gain PD controller as recommended by (Pastor et al., 2013; Schaal et al., 2007). Parameters of EDAs: ${\mathcal{L}}_{max}=2.5$ J.

For EDAs, both the controller with and without energy limitation were able to reach the goal (Figure 7(F) and (G)). However in the latter case, high accelerations of the end-effector occurred after the obstacle was removed (Figure 7(H)). As shown in Figure 8, λ is regulated to limit the total energy to be less than ${\mathcal{L}}_{max}$ . By including mechanical impedance as a separate primitive, the EDA approach was able to regulate the interactive dynamics between the robot and environment.

Figure 8.

Simulation results using Elementary Dynamic Actions (EDAs) with impedance modulation via energy regulation (equation (20)) (Figure 7) (Section 3.4.3). (Top) Time t versus λ of the controller. After the first contact, a sudden change in λ occurred. (Bottom) Time versus ${\mathcal{L}}_c$ of the controller. Since the virtual trajectory continued during contact, the decrease of λ limited ${\mathcal{L}}_c$ to a maximal value of 2.5 J. As expected, the total energy of the robot did not exceed ${\mathcal{L}}_{max}=2.5$ J.

3.5. Obstacle avoidance

We next consider obstacle avoidance while conducting a point-to-point discrete reaching movement presented in Section 3.3.3. We assume that the obstacle is fixed at location $\mathbf{o}\in {\mathbb{R}}^3$ ( ${\mathbb{R}}^2$ for the planar simulations), although the method can be easily generalized to nonstationary obstacles Andrews and Hogan (1983); Newman (1987).

3.5.1. Dynamic movement primitives

To avoid an obstacle located at $\mathbf{o}\in {\mathbb{R}}^3$ a coupling term C _t (t) is added to the transformation system (Hoffmann et al., 2009; Ijspeert et al., 2013; Saveriano et al., 2023): ${\tau }^2\ddot{\mathbf{p}}\left(t\right)=-{\alpha }_z{\beta }_z\left\{\mathbf{p}\left(t\right)-\mathbf{g}\right\}-{\alpha }_z\tau \stackrel{.}{\mathbf{p}}\left(t\right)+\mathbf{f}\left(s\left(t\right)\right)+{\mathbf{C}}_t\left(t\right)$

The coupling term C _t (t) is defined by: $\begin{array}{ll}{\mathbf{C}}_t\left(t\right)& =\gamma {\mathbf{R}}_c\left(t\right)\stackrel{.}{\mathbf{p}}\left(t\right)\theta \left(t\right)\exp \left(-\beta \theta \left(t\right)\right)\\ \theta \left(t\right)& =\arccos \left(\frac{{\left\{\mathbf{o}-\mathbf{p}\left(t\right)\right\}}^{\mathrm{T}}\stackrel{.}{\mathbf{p}}\left(t\right)}{\Vert \mathbf{o}-\mathbf{p}\left(t\right)\Vert \Vert \stackrel{.}{\mathbf{p}}\left(t\right)\Vert }\right)\end{array}$ (21)In these equations, γ and β are positive constants which correspond to the amplitude of the coupling term and its exponential decay rate, respectively; θ(t) is the angle between $\stackrel{.}{\mathbf{p}}\left(t\right)$ and o − p(t); $\Vert ·\Vert :{\mathbb{R}}^3\to {\mathbb{R}}_{\ge 0}$ is the Euclidean norm operator; R _c (t) ∈ SO(3) is a rotation matrix representing the ±90° rotation about axis $\mathbf{r}\left(t\right)\equiv \left\{\mathbf{o}-\mathbf{p}\left(t\right)\right\}\times \stackrel{.}{\mathbf{p}}\left(t\right)$ .

For spatial tasks, R _c (t) can be calculated using the Rodrigues’ formula (Lynch and Park, 2017; Murray et al., 1994): $\begin{array}{ll}{\mathbf{R}}_c\left(t\right)& ={\mathbf{I}}_3+\sin \left(\pm \frac{\pi }{2}\right)\left[\widehat{\mathbf{r}}\left(t\right)\right]+\left\{1-\hspace{0.17em}\cos \left(\pm \frac{\pi }{2}\right)\right\}{\left[\widehat{\mathbf{r}}\left(t\right)\right]}^2\\ & ={\mathbf{I}}_3\pm \left[\widehat{\mathbf{r}}\left(t\right)\right]+{\left[\widehat{\mathbf{r}}\left(t\right)\right]}^2\end{array}$ In this equation, $\widehat{\mathbf{r}}\left(t\right)$ is the normalization of r(t); $\left[\widehat{\mathbf{r}}\left(t\right)\right]$ is a skew-symmetric matrix form of $\widehat{\mathbf{r}}\left(t\right)$ (Murray et al., 1994). The + and − signs represent +90° and −90° rotation about axis $\widehat{\mathbf{r}}\left(t\right)$ , respectively. Intuitively, the coupling term forces the movement to move away from the obstacle (Ijspeert et al., 2013). For a planar task, R _c (t) is a skew-symmetric matrix with ±1 and ∓1 as off-diagonal terms.

3.5.2. Elementary dynamic actions

For EDAs, the idea resembles the method of obstacle avoidance using potential fields (Andrews and Hogan, 1983; Hjorth et al., 2020; Hogan, 1985; Khatib, 1985; Koditschek, 1987; Newman, 1987). A mechanical impedance which produces a repulsive force from the obstacle (i.e., a mechanical impedance with a point “repeller” at the obstacle location) is superimposed on the task-space impedance controller: $\begin{array}{l}{\mathbf{Z}}_{p,1}\left(t\right)={\mathbf{K}}_p\{{\mathbf{p}}_0\left(t\right)-\mathbf{p}\left(t\right))+{\mathbf{B}}_p\{{\stackrel{.}{\mathbf{p}}}_0\left(t\right)-\stackrel{.}{\mathbf{p}}\left(t\right)\}\\ {\mathbf{Z}}_{p,2}\left(t\right)=-\frac{k}{{\Vert \mathbf{o}-\mathbf{p}\left(t\right)\Vert }^n}\left(\mathbf{o}-\mathbf{p}\left(t\right)\right)\\ {\mathit{\tau }}_{in}\left(t\right)={\mathbf{J}}_p^{\mathrm{T}}\left(\mathbf{q}\right)\left\{{\mathbf{Z}}_{p,1}\left(t\right)+{\mathbf{Z}}_{p,2}\left(t\right)\right\}\end{array}$ (22)In these equations, k is a positive constant that determines the amplitude of the repulsive force; n is a positive integer. Intuitively, the impedance term Z_p,2(t) produces a “potential barrier” which repels the robot from the obstacle.

3.5.3. Simulation example

As in Section 3.3.3, we used the 2-DOF planar robot model to generate a goal-directed discrete movement in task-space coordinates. However, a stationary obstacle was located at o, blocking the path. The code script for this simulation is main_obstacle_avoidance.py.

As shown in Figure 9, both approaches successfully achieved the task. Nevertheless, differences between the two approaches were observed.

Figure 9.

A goal-directed discrete movement with obstacle avoidance for Dynamic Movement Primitives (DMPs, blue) and Elementary Dynamic Actions (EDAs, orange) (Section 3.5.3). (Left) Trajectory of the end-effector position. Initial end-effector position p _i and goal location g are depicted as black circle markers. (Right) Time-lapse of the movement for DMP (Top, blue) and EDA (Bottom, orange). Dotted lines are a minimum-jerk trajectory. Obstacle location o = [0.0, 1.14]m. Parameters of DMPs: γ = 300, β = −3 (equation (21)). Parameters of EDAs: k = 0.1 N/m⁵, n = 6 (equation (22)). Other parameters are identical to those in Figure 5.

DMPs considered the problem from a position control perspective. The coupling term C _t (t) (equation (21)) controlled the acceleration of the end-effector position p(t) for obstacle avoidance. As with Section 3.4, the presented example also highlighted the modular property of DMPs. Without modification, but simply reusing the learned weights from Section 3.3.3, DMPs superimposed coupling term C _t (t) onto the transformation systems to modify the learned trajectory. This example shows the property of online trajectory modulation of DMPs, which is an advantage over spline methods (Section 2.1.4).

In comparison, EDAs achieved obstacle avoidance without explicit path planning (Hogan, 1985). Instead, both goal-reaching and obstacle avoidance tasks were achieved by using the superposition principle of mechanical impedances (equation (10)). The modular property of the superposition principle enabled EDAs to divide the task into multiple sub-tasks, allocate an appropriate mechanical impedance for each sub-task, and then combine the mechanical impedances to solve the original task. For this case, mechanical impedance Z_p,1 was used for the goal-directed discrete movement planned in task-space coordinates and mechanical impedance Z_p,2 was used for obstacle avoidance. As with the weights of DMPs, Z_p,1 was simply reused from Section 3.3.3 without any modification.

While the modular property of EDAs simplified the approach, care is required for implementation. The method based on EDA is reminiscent of classic potential field methods (Khatib, 1986; Newman, 1987). Obstacle avoidance is achieved by superimposing a virtual (or artificial) elastic potential field that produces a repulsive force field emanating from the obstacle. Hence, the approach may encounter the known limitations of classic potential field methods. For instance, the end-effector of the robot may stall at a local energy minimum (Newman, 1987). To avoid this, a slight offset (2 cm) in the positive X direction of Figure 9 was used to “push” the robot toward the counterclockwise direction (Andrews and Hogan, 1983). Moreover, an inappropriate choice of controller parameters, such as an abrupt change of the virtual trajectory p₀(t), or relatively small damping matrices may cause undesired oscillation in the vicinity of the obstacle (Khatib, 1986; Muñoz Osorio et al., 2018). To address these problems, several methods, such as using time-varying potential fields, have been demonstrated. For brevity, further details are not reviewed here and the reader is referred to Andrews and Hogan (1983); Khatib (1986); Newman (1987); Muñoz Osorio et al. (2018); Hjorth et al. (2020); Becker et al. (2021); Huber et al. (2023).

3.6. Rhythmic movement

We next consider a method to generate a rhythmic, repetitive movement. For this example, we considered movements planned in both joint-space and task-space position.

3.6.1. Dynamic movement primitives

For DMPs, we used a canonical system and nonlinear forcing terms of a rhythmic movement (Section 2.1). From the generated q(t) (respectively p(t)), the process discussed in Section 3.2 (respectively Section 3.3) was used.

Compared to discrete movements (Section 3.2, 3.3), Imitation Learning is necessary to generate rhythmic movements, as f(s(t)) = 0 can only generate a goal-directed discrete movement (Section 2.1.4).

3.6.2. Elementary dynamic actions

For EDAs, we defined q₀(t) (respectively p₀(t)) to be a rhythmic movement which we aimed to follow. With this virtual trajectory, either a joint-space (equation (13)) or task-space impedance controller (equation (16)) was used.

3.6.3. Simulation example

We used the 2-DOF planar robot model from Section 3.2 and Section 3.3. The code scripts for the simulations are main_joint_rhythmic.py for joint-space and main_task_rhythmic.py for task-space.

The rhythmic movement in joint-space followed a sinusoidal trajectory. For DMPs and EDAs, both q _des (t) and q₀(t) were defined by: ${\mathbf{q}}_{des}\left(t\right),{\mathbf{q}}_0\left(t\right)={\mathbf{q}}_i+{\mathbf{q}}_A\hspace{0.17em}\sin \left({\omega }_{0}t\right)$ (23)In this equation, ${\mathbf{q}}_A\in {\mathbb{R}}^2$ is the amplitude of the rhythmic movement.

The rhythmic movement in task-space followed a circular trajectory. For DMPs and EDAs, both p _des (t) and p₀(t) were defined by: ${\mathbf{p}}_{des}\left(t\right),{\mathbf{p}}_0\left(t\right)=\mathbf{c}+\left[r_0\hspace{0.17em}\cos \left({\omega }_{0}t\right),r_0\hspace{0.17em}\sin \left({\omega }_{0}t\right)\right]$ (24)In this equation, $\mathbf{c}\in {\mathbb{R}}^2$ is the center location of the circle; r₀ is the radius of the circle.

As shown in Figure 10, both DMPs and EDAs successfully generated rhythmic movement in joint-space (Figures 10(A)–(D)) and task-space (Figures 10(E) and (F)). As discussed in Section 3.2.3 and Section 3.3.3, for DMPs, in principle, perfect tracking can be achieved in both joint-space and task-space. Given the period of the rhythmic movement (Section 2.1.4), rhythmic trajectories with arbitrary complexity can be learned using Imitation Learning. For EDAs, tracking error existed in both joint-space and task-space. Nevertheless, EDAs generated a rhythmic, repetitive movement without an inverse dynamics model and without solving the inverse kinematics. Moreover, as discussed in Section 3.2.3 and Section 3.3.3, the parameters of mechanical impedances can be chosen to reduce the tracking error (Figures 4 and 6).

Figure 10.

Rhythmic movements in (A, B, C, D) joint-space and (E, F) task-space for Dynamic Movement Primitives (DMPs, blue) and Elementary Dynamic Actions (EDAs, orange) (Section 3.6.3). (A, B) Rhythmic movement in joint-space and its (C, D) joint trajectories. (C, D) The black dashed lines which are perfectly overlapped with DMPs (blue lines) represent a sinusoidal trajectory with parameters: q _i = [0.5, 0.5] rad, q _A = [0.1, 0.3] rad, ω₀ = π rad/s (equation (23)). Parameters of rhythmic DMP: α _z = 10, β _z = 2.5, τ = 1.0/π, r = 1.0, N = 40, p = 100, c _i = 2π(i − 1)/N, h _i = N for i ∈ [1, 2, …, N]. Parameters of EDA are identical to those in Figure 3. (E, F) Rhythmic movement in task-space. The black dashed lines represent a circular trajectory with parameters: c = [0, 1.4142]m, r₀ = 0.5 m, ω₀ = π rad/s (equation (24)). Parameters of EDAs: K _p = 90I₂ N/m, B _p = 60I₂ N⋅s/m. Higher stiffness and damping values than Figure 5 were chosen to reduce the tracking error. Parameters of rhythmic DMP are identical to those in (A, B). Consistent with Figures 3 and 5, for DMP, perfect tracking was achieved while for EDA, a non-negligible tracking error existed.

3.7. Combination of discrete and rhythmic movements

We next designed a controller to generate a combination of discrete and rhythmic movements (Hogan and Sternad, 2007). For this example, we considered a movement planned in both joint-space and task-space positions.

3.7.1. Dynamic movement primitives

For DMPs, the canonical system and nonlinear forcing term are different for discrete and rhythmic movements (Section 2.1.1, 2.1.2). Hence, both rhythmic and discrete DMPs cannot be directly combined. Instead, the discrete DMP generates a time-changing goal g(t) for the rhythmic DMP (Degallier et al., 2006, 2007, 2008): $\begin{array}{l}{\tau }_1^2\ddot{\mathbf{q}}\left(t\right)=-{\alpha }_{z,1}{\beta }_{z,1}\left\{\mathbf{q}\left(t\right)-\mathbf{g}\left(t\right)\right\}-{\alpha }_{z,1}{\tau }_1\stackrel{.}{\mathbf{q}}\left(t\right)+\mathbf{f}\left(s\left(t\right)\right)\\ s\left(t\right)=\frac{t}{{\tau }_1}\mathrm{mod}\hspace{0.17em}2\pi \\ {\tau }_2^2\ddot{\mathbf{g}}\left(t\right)=-{\alpha }_{z,2}{\beta }_{z,2}\left\{\mathbf{g}\left(t\right)-{\mathbf{g}}_0\right\}-{\alpha }_{z,2}{\tau }_2\stackrel{.}{\mathbf{g}}\left(t\right)\end{array}$ (25)

The first two equations represent rhythmic DMPs, and the last equation represents discrete DMPs without a nonlinear forcing term and a canonical system. g₀ is the discontinuous change of the goal location to which goal g(t) converges. Note that for control in task-space position, q(t) terms in the equation are substituted for p(t). From the generated q(t) (respectively p(t)), the process discussed in Section 3.2 (respectively Section 3.3) was used.

3.7.2. Elementary dynamic actions

For EDAs, we simply add submovements q_0,sub(t) and oscillations q_0,osc(t) to define the virtual trajectory q₀(t): ${\mathbf{q}}_0\left(t\right)={\mathbf{q}}_{0,sub}\left(t\right)+{\mathbf{q}}_{0,osc}\left(t\right)$

Note that for control in task-space, q terms in the equation are substituted with p. With this virtual trajectory, either a joint-space (equation (13)) or a task-space impedance controller (equation (16)) is used.

3.7.3. Simulation example

Using the 2-DOF robot model in Section 3.2, the goal was to generate a combination of discrete and rhythmic movements both in joint-space and task-space. The code scripts for the simulations are main_joint_discrete_and_rhythmic.py for joint-space and main_task_discrete_and_rhythmic.py for task-space position.

For q _des (t) (respectively p _des (t)) of DMPs, the sinusoidal trajectory, equation (23) (respectively circular trajectory, equation (24)) was represented by rhythmic DMPs, and g(t) of q _des (t) (respectively p _des (t)) was represented by discrete DMPs with g₀ discretely changing from q _i to q _f (respectively p _i to p _f ). For control in joint-space, g₀ changed between q _i = [0.5, 0.5]rad and q _f = [1.5, 1.5]rad (Figures 11(A)–11(F)). For control in the task-space position, g₀ changed between p _i = [−0.47, 0.9]m and p _f = [0.53, 0.9]m (Figures 11(G)–11(J)). The time at which the i-th discrete movement started was 5i-2 seconds for $i\in \{1,2\}$

Figure 11.

A combination of discrete and rhythmic movements in (A–F) joint-space and (G–J) task-space for Dynamic Movement Primitives (DMPs, blue) and Elementary Dynamic Actions (EDAs, orange) (Section 3.7.3). (C, E, I) Black filled lines represent the discrete change of g₀ of DMP. Black dotted lines represent g(t), which follows a response of a (critically damped) second-order linear system. (D, F, J) Black dotted lines represent q₀(t) and p₀(t) for joint-space and task-space, respectively. (A–F) Parameters of the discrete and rhythmic movements in joint-space: ω₀ = π rad/s, q _A = [0.1, 0.3] rad, q _i = [0.5, 0.5] rad, q _f = [1.5, 1.5] rad. (G–J) Parameters of the discrete and rhythmic movements in task-space: ω₀ = π rad/s, r₀ = 0.3 m, p _i = [−0.47, 0.9]m, p _f = [0.53, 0.9]m. In joint-space (respectively task-space), multiple discrete movements in opposite directions are generated by switching the values of q _i (respectively p _i ) and q _f (respectively p _f ), with time offset t _off = 5i-2.0, for the i-th discrete movement, $i\in \{1,2\}$ . Parameters of discrete DMPs in both joint-space and task-space: τ₂ = 1.0, α_z,2 = 10, β_z,2 = 2.5. Other parameters of DMP are identical to those in Figure 10. Parameters of EDA are identical to those in Figure 10. For EDA, given a single impedance operator, discrete and rhythmic movements can be directly combined at the level of the virtual trajectory.

For q₀(t) of EDAs, a minimum-jerk trajectory (equation (14)) was combined with a sinusoidal trajectory (equation (23)): $\begin{array}{l}{\mathbf{q}}_0\left(t\right)={\mathbf{q}}_{0,sub}\left(t\right)+{\mathbf{q}}_{0,osc}\left(t\right)\\ {\mathbf{q}}_{0,sub}\left(t\right)=\{\begin{array}{ll}\mathbf{0}& 0\le t<t_{off}\\ {\mathbf{q}}_i+\left({\mathbf{q}}_f-{\mathbf{q}}_i\right)f_{MJT}\left(t\right)& t_{off}\le t<T+t_{off}\\ {\mathbf{q}}_f& T+t_{off}\le t\end{array}\\ {\mathbf{q}}_{0,osc}\left(t\right)={\mathbf{q}}_A\hspace{0.17em}\sin \left({\omega }_{0}t\right)\\ f_{MJT}\left(t\right)=10{\left(\frac{t-t_{off}}{T}\right)}^3-15{\left(\frac{t-t_{off}}{T}\right)}^4+6{\left(\frac{t-t_{off}}{T}\right)}^5\end{array}$ In this equation, t _off > 0 is a time offset for the submovement. For q _des (t) of DMPs, the sinusoidal trajectory was represented by rhythmic DMPs, and g(t) of q _des (t) was represented by discrete DMPs with g₀ discretely changing from q _i to q _f . For submovement q_0,sub(t), the parameters were x_0,i = q _i , v = q _f − q _i , and f _σ (t) = f _MJT (t) for t _off ≤ t ≤ T + t _off (equation (8)).

For p₀(t) of EDAs, a minimum-jerk trajectory (equation (14)) was combined with a circular trajectory (equation (24)): $\begin{array}{l}{\mathbf{p}}_0\left(t\right)={\mathbf{p}}_{0,sub}\left(t\right)+{\mathbf{p}}_{0,osc}\left(t\right)\\ {\mathbf{p}}_{0,sub}\left(t\right)=\{\begin{array}{ll}\mathbf{0}& 0\le t<t_{off}\\ {\mathbf{p}}_i+\left({\mathbf{p}}_f-{\mathbf{p}}_i\right)f_{MJT}\left(t\right)& t_{off}\le t<T+t_{off}\\ \mathbf{g}& T+t_{off}\le t\end{array}\\ {\mathbf{p}}_{0,osc}\left(t\right)=\left[r_0\hspace{0.17em}\cos \left({\omega }_{0}t\right),r_0\hspace{0.17em}\sin \left({\omega }_{0}t\right)\right]\end{array}$

For submovement p_0,sub(t), the parameters were x_0,i = p _i , v = p _f − p _i , and f _σ (t) = f _MJT (t) for t _off ≤ t ≤ T + t _off (equation (8)).

As shown in Figure 11, both approaches successfully produced a combination of discrete and rhythmic movements in joint-space and task-space. However, since DMPs separate the canonical system and nonlinear forcing terms for discrete and rhythmic movements (Section 2.1.1, 2.1.2), merging the two movements was not straightforward. DMPs circumvented this issue by assigning the time-changing goal g(t) of rhythmic DMPs, with discrete DMPs without the nonlinear forcing terms input.

On the other hand, for EDAs, given a single impedance operator, discrete and rhythmic movements were directly combined at the level of the virtual trajectory (i.e., the forward path dynamics) (Figure 2). With modest parameter tuning, the discrete and rhythmic movements used in Section 3.2.3 and Section 3.6.3 were reused and combined. This approach intuitively provides convenience in practical implementation and also emphasizes the modularity of EDAs at a kinematic level.

3.8. Sequence of discrete movements

We next consider designing a controller to generate a sequence of discrete movements planned in task-space position. For this, we show how the controller generates a movement in response to a sudden change of goal location.

3.8.1. Dynamic movement primitives

With the controller introduced in Section 3.3.1, an additional differential equation for the time-varying goal location g(t) was added (Ijspeert et al., 2013): $\begin{array}{l}{\tau }^2\ddot{\mathbf{p}}\left(t\right)=-{\alpha }_z{\beta }_z\left\{\mathbf{p}\left(t\right)-\mathbf{g}\left(t\right)\right\}-{\alpha }_z\tau \stackrel{.}{\mathbf{p}}\left(t\right)+\mathbf{f}\left(s\left(t\right)\right)\\ \tau \stackrel{.}{\mathbf{g}}\left(t\right)={\alpha }_g\left\{{\mathbf{g}}_0-\mathbf{g}\left(t\right)\right\}\end{array}$ (26)In these equations, α _g is a positive constant; g₀ is a discontinuous change of the goal location. In other words, the goal location g is first-order low-pass filtered. This is used to avoid a discontinuous jump in the acceleration $\ddot{\mathbf{p}}\left(t\right)$ of the transformation system.

Note that both p(t) and g(t) comprise a third-order linear system with a nonlinear forcing term and α _g g₀/τ as inputs. While any positive values of α _g can be used, for τ = 1 and β _z = α _z /4, α _g is often chosen to be α _g = α _z /2 such that the third-order linear system has three repeated eigenvalues (Nemec and Ude, 2012; Schaal et al., 2005).

Note that g(t) has a closed-form solution. Without loss of generality, if the goal location changes from g₀ = g _old to g₀ = g _new at time t = 0, with initial condition g(0) = g _old : $\mathbf{g}\left(t\right)={\mathbf{g}}_{new}+\left({\mathbf{g}}_{old}-{\mathbf{g}}_{new}\right)\exp \left(-\frac{{\alpha }_g}{\tau }t\right)$

Hence, a sequence of finite submovements can be generated by discrete changes of the goal location from g₀ = g _old to g₀ = g _new .

While the first discrete movement can follow an arbitrary trajectory using Imitation Learning, with this formulation (Nemec and Ude, 2012), the subsequent discrete movements cannot, but follow the motion of a stable third-order linear system converging to the corresponding g₀ value.

3.8.2. Elementary dynamic actions

With the first-order position task-space impedance controller (Equation (16)), a sequence of discrete movements is generated by sequencing multiple submovements (Flash and Henis, 1991): ${\stackrel{.}{\mathbf{p}}}_0\left(t\right)={\displaystyle \sum }_i{\mathbf{v}}_i{\widehat{\sigma }}_i\left(t-t_i\right)$

The amplitude of the ith basis function v _i is chosen to reach the goal of the ith submovement, starting from the goal of the previous submovement. Note that each submovement can use a different basis function $\widehat{\sigma }$ .

3.8.3. Simulation example

The simulation example in Figure 12 reproduced the experiment conducted in Flash and Henis (1991). Let the goal location of the first discrete movement be g _old . A new goal location, g _new suddenly appeared at time t _g . The task was to modulate the controller to eventually reach g _new . We assumed that the new goal position g _new was immediately measured at time t _g . The simulation example used the 2-DOF robot model of Section 3.2. The code script for the simulation is main_sequencing.py

Figure 12.

A sequence of discrete movements for (A, B, D, E) Dynamic Movement Primitives (DMPs, blue) and (A, C, D, E) Elementary Dynamic Actions (EDAs, orange) (Section 3.8.3). (A) Trajectory of the end-effector position. (B, C) Time-lapse of the movement of the 2-DOF robot model. (D, E) Time versus X- and Y-coordinates of the end-effector trajectory. (A, B, C) The first movement headed toward g _old (square grey marker), but at time t = t _g (D, E), the target switched location to g _new (square black marker) which necessitated a second movement. Parameters of DMPs (for the first discrete movement): α _g = 1.0, τ = 1.0. Other parameters are identical to those in Figure 5. Parameters of EDAs: T₁ = 1.0s, T₂ = 1.0s, t _g = 0.5s. g _old = [−0.7, 1.22]m, g _new = [0.8, 1.72]m, p _i = [0.0, 0.52]m. Other parameters are identical to those in Figure 5. (D, E) For the new goal location g _new , EDA achieved faster convergence than DMP. For DMP, using the formulation of Nemec and Ude (2012), the second movement followed a response of a third-order linear system. On the other hand for EDA, the basis function of the second movement can be freely chosen. Hence, faster convergence to g _new was achieved for EDA.

For DMPs, the first discrete movement followed a minimum-jerk trajectory with goal location g _old . The subsequent discrete movement was generated by discretely changing g₀ = g _old to g₀ = g _new at t = t _g : ${\mathbf{g}}_0\left(t\right)=\{\begin{array}{ll}{\mathbf{g}}_{old}& 0\le t<t_g\\ {\mathbf{g}}_{new}& t_g\le t\end{array}$

For EDAs, the minimum-jerk trajectory was used for the basis function of each submovement p₀(t): $\begin{array}{l}{\mathbf{p}}_0\left(t\right)={\mathbf{p}}_{0,1}\left(t\right)+{\mathbf{p}}_{0,2}\left(t\right)\\ {\mathbf{p}}_{0,1}\left(t\right)=\{\begin{array}{ll}{\mathbf{p}}_i+\left({\mathbf{g}}_{old}-{\mathbf{p}}_i\right)f_{MJT,1}\left(t\right)& 0\le t<T_1\\ {\mathbf{g}}_{old}& T_1\le t\end{array}\\ {\mathbf{p}}_{0,2}\left(t\right)=\{\begin{array}{ll}\mathbf{0}& 0\le t<t_g\\ \left({\mathbf{g}}_{new}-{\mathbf{g}}_{old}\right)f_{MJT,2}\left(t\right)& t_g\le t<t_g+T_2\\ {\mathbf{g}}_{new}-{\mathbf{g}}_{old}& t_g+T_2\le t\end{array}\\ f_{MJT,1}\left(t\right)=10{\left(\frac{t}{T_1}\right)}^3-15{\left(\frac{t}{T_1}\right)}^4+6{\left(\frac{t}{T_1}\right)}^5\\ f_{MJT,2}\left(t\right)=10{\left(\frac{t-t_g}{T_2}\right)}^3-15{\left(\frac{t-t_g}{T_2}\right)}^4+6{\left(\frac{t-t_g}{T_2}\right)}^5\end{array}$ In these equations, subscripts 1 and 2 denote the first and second submovements, respectively; For submovement p_0,1(t), the parameters were x_0,i = p _i , v = g _old − p _i , and f _σ (t) = f_MJT,1(t) for 0 ≤ t ≤ T₁; for submovement p_0,2(t), the parameters were x_0,i = 0, v = g _new − g _old , and f _σ (t) = f_MJT,2(t) for t _g ≤ t ≤ T₂ + t _g (equation (8)). With this p₀(t), the controller introduced in Section 3.3.2 was used.

As shown in Figure 12, for both approaches, p(t) converged to the new goal location g _new . Nevertheless, it should be noted that the modular property of EDAs at the kinematics level enabled a smooth integration of the second movement without any modification of the first one.

3.9. Managing kinematic redundancy

We next consider designing a controller to generate a goal-directed (or sequence of) discrete movement(s) of the end-effector’s position for a kinematically redundant robot. By definition, kinematic redundancy occurs when a Jacobian matrix has a null space (Siciliano, 1990). Hence, infinitely many joint velocity solutions exist to produce a desired end-effector velocity. While kinematic redundancy provides significant challenges, additional control objectives can be achieved by exploiting kinematic redundancy. Examples include obstacle avoidance (Baillieul, 1986; Maciejewski and Klein, 1985), joint limit avoidance (Hjorth et al., 2020; Liegeois, 1977), and minimization of instantaneous power during movement (Klein and Huang, 1983).

3.9.1. Dynamic movement primitives

For DMPs, a feedback controller is employed to manage kinematic redundancy. Multiple feedback control methods exist and can be divided into three categories: velocity-based control, acceleration-based control, and force-based control (Nakanishi et al., 2005, 2008). Among these methods, we used a “velocity-based control without joint-velocity integration” (Nakanishi et al., 2008; Pastor et al., 2009).

Let p _des (t) be the desired end-effector trajectory with duration T and p _des (T) = g. By generating p _des (t), ${\stackrel{.}{\mathbf{p}}}_{des}\left(t\right)$ , ${\ddot{\mathbf{p}}}_{des}\left(t\right)$ as shown in Section 3.3, a reference end-effector velocity ${\stackrel{.}{\mathbf{p}}}_r\left(t\right)$ and its corresponding reference joint velocity ${\stackrel{.}{\mathbf{q}}}_r\left(t\right)$ are defined: $\begin{array}{l}{\stackrel{.}{\mathbf{p}}}_r\left(t\right)={\stackrel{.}{\mathbf{p}}}_{des}\left(t\right)+{\mathbf{\Lambda }}_p\left({\mathbf{p}}_{des}\left(t\right)-\mathbf{p}\left(t\right)\right)\\ {\stackrel{.}{\mathbf{q}}}_r\left(t\right)={\mathbf{J}}_p^{†}\left(\mathbf{q}\right){\stackrel{.}{\mathbf{p}}}_r\left(t\right)\end{array}$ (27)In these equations, ${\mathbf{J}}_p^{†}\left(\mathbf{q}\right)$ denotes the Moore-Penrose pseudo inverse, which is defined by ${\mathbf{J}}_p^{†}\left(\mathbf{q}\right)={\mathbf{J}}_p^{\mathrm{T}}\left(\mathbf{q}\right){\left\{{\mathbf{J}}_p\left(\mathbf{q}\right){\mathbf{J}}_p^{\mathrm{T}}\left(\mathbf{q}\right)\right\}}^{-1}$ (Klein and Huang, 1983; Penrose, 1955); ${\mathbf{\Lambda }}_p\in {\mathbb{R}}^{3\times 3}$ ( ${\mathbb{R}}^{2\times 2}$ for planar task) is a symmetric positive definite matrix. Note that for ${\stackrel{.}{\mathbf{q}}}_r\left(t\right)$ , an additional term which projects an arbitrary joint velocity vector onto the null space of the Jacobian matrix J _p (q) can be added (Liegeois, 1977). However, this additional term was omitted for brevity.

Accordingly, the reference end-effector acceleration ${\ddot{\mathbf{p}}}_r\left(t\right)$ and its corresponding reference joint acceleration ${\ddot{\mathbf{q}}}_r\left(t\right)$ are defined: $\begin{array}{ll}{\ddot{\mathbf{p}}}_r\left(t\right)& ={\ddot{\mathbf{p}}}_{des}\left(t\right)+{\mathbf{\Lambda }}_p\left({\stackrel{.}{\mathbf{p}}}_{des}\left(t\right)-\stackrel{.}{\mathbf{p}}\left(t\right)\right)\\ {\ddot{\mathbf{q}}}_r\left(t\right)& ={\mathbf{J}}_p^{†}\left(\mathbf{q}\right)\left({\ddot{\mathbf{p}}}_r\left(t\right)-{\stackrel{.}{\mathbf{J}}}_p\left(\mathbf{q}\right)\stackrel{.}{\mathbf{q}}\left(t\right)\right)\\ & ={\mathbf{J}}_p^{†}\left(\mathbf{q}\right)\left\{{\ddot{\mathbf{p}}}_{des}\left(t\right)+{\mathbf{\Lambda }}_p\left[{\stackrel{.}{\mathbf{p}}}_{des}\left(t\right)-\stackrel{.}{\mathbf{p}}\left(t\right)\right]-{\stackrel{.}{\mathbf{J}}}_p\left(\mathbf{q}\right)\stackrel{.}{\mathbf{q}}\left(t\right)\right\}\end{array}$

From these values, τ _in (t) is defined by: ${\mathit{\tau }}_{in}\left(t\right)=\mathbf{M}\left(\mathbf{q}\right){\ddot{\mathbf{q}}}_r\left(t\right)+\mathbf{C}\left(\mathbf{q},\stackrel{.}{\mathbf{q}}\right){\stackrel{.}{\mathbf{q}}}_r\left(t\right)-{\mathbf{\Lambda }}_q\left\{\stackrel{.}{\mathbf{q}}\left(t\right)-{\stackrel{.}{\mathbf{q}}}_r\left(t\right)\right\}$ In this equation, ${\mathbf{\Lambda }}_q\in {\mathbb{R}}^{n\times n}$ is a symmetric positive definite matrix.

Note that this controller, suggested in Nakanishi et al. (2008), is equivalent to the sliding mode feedback controller introduced by Slotine and Li (1987). It was shown by Slotine and Li (1987) that p(t) asymptotically converges to p _des (t). Moreover, with this feedback controller, one can resolve kinematic singularity using a Damped Least-squares inverse (Section 3.3.3) (Chiaverini et al., 1994; Nakamura and Hanafusa, 1986).

3.9.2 Elementary dynamic actions

For EDAs, kinematic redundancy of the robot manipulator can be managed by superimposing multiple mechanical impedances (Hermus et al., 2021; Verdi, 2019). In detail, a first-order position task-space impedance controller (equation (16)) can be superimposed with a joint-space controller that implements joint damping (equation (13)): $\begin{array}{l}{\mathbf{Z}}_p\left(t\right)={\mathbf{K}}_p\left\{{\mathbf{p}}_0\left(t\right)-\mathbf{p}\left(t\right)\right\}+{\mathbf{B}}_p\left\{{\stackrel{.}{\mathbf{p}}}_0\left(t\right)-\stackrel{.}{\mathbf{p}}\left(t\right)\right\}\\ {\mathbf{Z}}_q\left(t\right)=-{\mathbf{B}}_q\stackrel{.}{\mathbf{q}}\left(t\right)\\ {\mathit{\tau }}_{in}\left(t\right)={\mathbf{J}}_p^{\mathrm{T}}\left(\mathbf{q}\right){\mathbf{Z}}_p\left(t\right)+{\mathbf{Z}}_q\left(t\right)\end{array}$ (28)

As in Section 3.2.2 and Section 3.3.2, with constant symmetric positive definite matrices of K _p , B _p , B _q , a goal directed discrete movement can be generated by setting p₀(t) to be a submovement which ends at goal location g.

The stability of this controller is shown in Arimoto et al. (2005a, 2005b); Arimoto and Sekimoto (2006); Lachner (2022), where an asymptotic convergence of $\stackrel{.}{\mathbf{q}}\left(t\right)\to \mathbf{0}$ and p(t) → g is guaranteed for discrete movement p₀(t). This also implies an asymptotic convergence of q(t) to one of the infinite number of solutions that satisfies g = h _p (q(t)).

Like the controller in equation (16), this controller does not involve an inversion of the Jacobian matrix. Moreover, explicitly solving the inverse kinematics and an inverse dynamics model are not required. Hence, the approach remains stable near kinematic singularities.

For the desired discrete movement, we used a damping joint-space impedance operator to reduce the joint motions in the nullspace of J _p (q) (equation (28)). Note that for repetitive rhythmic movements in task-space (Klein and Huang, 1983), this controller might result in non-negligible drift in joint space (Mussa-Ivaldi and Hogan, 1991). If this resulted in control problems, for example, violation of joint limits, one can augment the damping impedance operator with a stiffness term K _q (q₀(t) − q(t)), which will eliminate the joint-space drift and enable a stable equilibrium configuration. For brevity, this was not considered here.

Superimposing joint-space and task-space impedances can yield task conflicts (Schettino et al., 2021), unless the virtual joint configuration to which the joint stiffness is connected is defined at the desired goal location (Hermus et al., 2021). Often, the task conflict is resolved by using null-space projection methods, as suggested by Khatib (1995). However, it is important to note that the resultant controller violates passivity (Lachner, 2022). Alternatively, it has been shown that, with sufficiently large null-space dimension, the task conflict is minimized or even eliminated (Hermus et al., 2021).

3.9.3. Simulation example

Consider a 5-DOF planar serial-link robot model, where each link consists of a single uniform slender bar with mass and length of 1 kg and 1m, respectively. With this robot model, we generated a single (or sequence of) discrete movement(s). As in Section 3.3.3, a minimum-jerk trajectory was used. For the sequence of discrete movements, the trajectories of Section 3.8.3 were used. The code scripts for the simulations are main_redundant_discrete.py for single discrete movement and main_redundant_sequencing.py for sequence of discrete movements.

As shown in Figure 13, both approaches were able to achieve goal-directed discrete movements. For DMPs, the approach used a feedback controller with reference trajectories generated by DMPs to manage kinematic redundancy (Nakanishi et al., 2008; Slotine and Li, 1991). For EDAs, the controller simply reused the task-space impedance controller in Section 3.3.2 and combined it with the impedance controller (with K _q = 0) of Section 3.2.2. As demonstrated in Section 3.5, this example again emphasizes the modular property of EDAs.

Figure 13.

Managing kinematic redundancy using Dynamic Movement Primitives (DMPs, blue) and Elementary Dynamic Actions (EDAs, orange) (Section 3.9.3). (A, B, C) A single discrete movement. (D, E, F) A sequence of discrete movements. The first movement headed toward g _old (square grey marker), but at time t = t _g , target switched location to g _new (square black marker) which necessitated a second movement. (A, B) The dashed black lines are a minimum-jerk trajectory. (C, F) End-effector trajectories. (A, B, C) Parameters of the minimum-jerk trajectory: p _i = [0, 3]m, g = [3, 3]m, T = 2.0s. (D, E, F) Parameters of the minimum-jerk trajectories: p _i = [−1.62, 0.76]m, g _old = [−3.62, 1.76]m, g _new = [2.38, 3.26]m, T₁ = 2.0s, T₂ = 3.0s, t _g = 1.0. Parameters of DMPs: Λ _p = 80I₂, Λ _q = 100I₅, τ = T = T₁. Other parameters are identical to those in Figure 3. Parameters of EDAs: K _p = 300I₂ N/m, B _p = 100I₂ N⋅s/m, B _q = 30I₅ N⋅m⋅s/rad.

3.10. Managing kinematic redundancy while controlling position and orientation

We next consider controlling both the position and orientation of the robot’s end-effector, while managing kinematic redundancy. As with Section 3.9, we consider a discrete movement in task-space.

3.10.1. Dynamic movement primitives

Various DMPs have been proposed to represent task-space orientation (Abu-Dakka et al., 2015; Pastor et al., 2011; Saveriano et al., 2019, 2023; Ude et al., 2014). Koutras and Doulgeri (2020) suggested a formulation using unit quaternions, which addressed the limitations of prior methods.

In detail, DMPs to represent task-space orientation were defined by (Koutras and Doulgeri, 2020): $\begin{array}{l}{\tau }^2\stackrel{..}{\mathbf{e}}\left(t\right)=-{\alpha }_z{\beta }_z\mathbf{e}\left(t\right)-{\alpha }_z\tau \stackrel{.}{\mathbf{e}}\left(t\right)+\mathbf{f}\left(s\left(t\right)\right)\\ \tau \dot{s}\left(t\right)=-{\alpha }_{s}s\left(t\right)\end{array}$ (29)In these equations, $\mathbf{e}\left(t\right)\equiv 2\mathrm{Log}({\overrightarrow{\mathbf{q}}}_g\otimes \overrightarrow{\mathbf{q}}\left(t\right)*)\in {\mathrm{ℝ}}^3$ denotes the orientation error vector; $\overrightarrow{\mathbf{q}}\left(t\right),{\overrightarrow{\mathbf{q}}}_g\in {\mathbb{H}}_1$ denotes the current and goal orientations of the robot’s end-effector; $\mathrm{Log}:{\mathbb{H}}_1\to {\mathbb{R}}^3$ is a Logarithm Map (Koutras and Doulgeri, 2020; Kuipers, 1999); $\otimes :{\mathbb{H}}_1\times {\mathbb{H}}_1\to {\mathbb{H}}_1$ denotes the quaternion multiplication operator; $*:{\mathbb{H}}_1\to {\mathbb{H}}_1$ denotes the quaternion conjugate operator.

In contrast to joint-space (equation (12)) and task-space position (equation (15)) trajectories, the following nonlinear forcing term was used for task-space orientation (equation (29)) (Koutras and Doulgeri, 2020): $\mathbf{f}\left(s\left(t\right)\right)=2\text{diag}\left(\mathrm{Log}\left({\overrightarrow{\mathbf{q}}}_g\otimes {\overrightarrow{\mathbf{q}}}_i^{*}\right)\right)\frac{{\sum }_{i=1}^N{\mathbf{w}}_i{\varphi }_i\left(s\left(t\right)\right)}{{\sum }_{i=1}^N{\varphi }_i\left(s\left(t\right)\right)}s\left(t\right)$ In this equation, $\text{diag}(·):{\mathbb{R}}^3\to {\mathbb{R}}^{3\times 3}$ defines a diagonal matrix with its elements defined by the arguments; ${\overrightarrow{\mathbf{q}}}_i\in {\mathbb{H}}_1$ denotes the initial orientation of the robot’s end-effector; ${\mathbf{w}}_i\in {\mathbb{R}}^3$ is the weight array of the ith basis function. Compared to prior examples, the nonlinear forcing term f(s(t)) instead includes the term $2\text{diag}\left(\mathrm{Log}\left({\overrightarrow{\mathbf{q}}}_g\otimes {\overrightarrow{\mathbf{q}}}_i^{*}\right)\right)$ rather than g − q _i for joint-space trajectory (respectively g − p _i for task-space position trajectory).

The formulation of this transformation system (equation (29)) is identical to a three-dimensional DMP. Hence, given P samples of ${\ddot{\mathbf{e}}}_{des}\left(t\right),{\stackrel{.}{\mathbf{e}}}_{des}\left(t\right),{\mathbf{e}}_{des}\left(t\right)$ calculated from ${\overrightarrow{\mathbf{q}}}_{des}\left(t\right)$ and its higher derivatives, one can use Imitation Learning (or Learning from Demonstration) using Locally Weighted regression (or Linear Least-Square regression) (Section 2.1.4) (Koutras and Doulgeri, 2020).

Once e _des (t) is generated by Imitation Learning, the desired orientation of the end-effector, ${\overrightarrow{\mathbf{q}}}_{des}\left(t\right)$ is reconstructed from: ${\overrightarrow{\mathbf{q}}}_{des}\left(t\right)=\mathrm{Exp}{\left(\frac{1}{2}{\overrightarrow{\mathbf{e}}}_{des}\left(t\right)\right)}^{*}\otimes {\overrightarrow{\mathbf{q}}}_g$ (30)In this equation, $\mathrm{Exp}:{\mathbb{R}}^3\to {\mathbb{H}}_1$ is the Exponential Map (Koutras and Doulgeri, 2020).

Given ${\overrightarrow{\mathbf{q}}}_{des}\left(t\right)$ , and p _des (t) generated by separate DMPs (Section 3.3.1, equation (15)), tracking control of both task-space position and orientation can be accomplished using the following feedback controller (Ceccarelli, 2008): $\begin{array}{l}{\mathit{\tau }}_{in}\left(t\right)=\mathbf{M}\left(\mathbf{q}\right)\left({\stackrel{.}{\mathbf{x}}}_{des}\left(t\right)+{\mathbf{\Lambda }}_q\stackrel{\sim }{\mathbf{x}}\left(t\right)\right)+\mathbf{C}\left(\mathbf{q},\stackrel{.}{\mathbf{q}}\right)\stackrel{.}{\mathbf{q}}\left(t\right)\\ \stackrel{\sim }{\mathbf{x}}\left(t\right)={\mathbf{x}}_{des}\left(t\right)-\stackrel{.}{\mathbf{q}}\left(t\right)\\ {\mathbf{x}}_{des}\left(t\right)={\mathbf{J}}^{†}\left(\mathbf{q}\right)\left[\begin{array}{c}{\stackrel{.}{\mathbf{p}}}_{des}\left(t\right)+{\mathbf{\Lambda }}_p\left({\mathbf{p}}_{des}\left(t\right)-\mathbf{p}\left(t\right)\right)\\ {\mathbf{\omega }}_{des}\left(t\right)+{\mathbf{\Lambda }}_r\stackrel{\sim }{ϵ}\left(t\right)\end{array}\right]\\ \stackrel{\sim }{ϵ}\left(t\right)=\eta \left(t\right){\mathbf{ϵ}}_{des}\left(t\right)-{\eta }_{des}\left(t\right)\mathbf{ϵ}\left(t\right)+\left[\mathbf{ϵ}\left(t\right)\right]{\mathbf{ϵ}}_{des}\left(t\right)\end{array}$ (31)In these equations, ${\mathbf{\Lambda }}_r\in {\mathbb{R}}^{3\times 3}$ and ${\mathbf{\Lambda }}_q\in {\mathbb{R}}^{n\times n}$ are symmetric positive definite matrices; $\overrightarrow{\mathbf{q}}\left(t\right)=\left(\eta \left(t\right),\mathbf{ϵ}\left(t\right)\right)$ and ${\overrightarrow{\mathbf{q}}}_{des}\left(t\right)=\left({\eta }_{des}\left(t\right),{\mathbf{ϵ}}_{des}\left(t\right)\right)$ ; $\mathbf{J}\left(\mathbf{q}\right)=\left[{\mathbf{J}}_p\left(\mathbf{q}\right),{\mathbf{J}}_r\left(\mathbf{q}\right)\right]\in {\mathbb{R}}^{6\times n}$ ; ${\mathbf{\omega }}_{des}\left(t\right)\in {\mathbb{R}}^3$ is the desired angular velocity of the end-effector, which can be derived from ${\overrightarrow{\mathbf{q}}}_{des}\left(t\right)$ and its higher derivatives (Koutras and Doulgeri, 2020).

While an analytical form of ${\stackrel{.}{\mathbf{x}}}_{des}\left(t\right)$ can be derived (Ceccarelli, 2008), for simplicity, numerical differentiation may be used.

3.10.2. Elementary dynamic actions

Using EDAs, the superposition principle of mechanical impedances (equation (10)) may be exploited by simply adding an orientation task-space impedance operator onto equation (28). In detail, the following impedance operator for orientation, Z _r (Hermus et al., 2021) can be superimposed: $\begin{array}{l}{\mathbf{Z}}_r\left(t\right)=k_r\mathbf{R}\left(t\right)\mathrm{Log}\left({\mathbf{R}}^{\mathrm{T}}\left(t\right){\mathbf{R}}_0\left(t\right)\right)-b_r\mathit{\omega }\left(t\right)\\ {\mathbf{Z}}_p\left(t\right)={\mathbf{K}}_p\left\{{\mathbf{p}}_0\left(t\right)-\mathbf{p}\left(t\right)\right\}+{\mathbf{B}}_p\left\{{\stackrel{.}{\mathbf{p}}}_0\left(t\right)-\stackrel{.}{\mathbf{p}}\left(t\right)\right\}\\ {\mathbf{Z}}_q\left(t\right)=-{\mathbf{B}}_q\stackrel{.}{\mathbf{q}}\left(t\right)\\ {\mathit{\tau }}_{in}\left(t\right)={\mathbf{J}}_r^{\mathrm{T}}\left(\mathbf{q}\right){\mathbf{Z}}_r\left(t\right)+{\mathbf{J}}_p^{\mathrm{T}}\left(\mathbf{q}\right){\mathbf{Z}}_p\left(t\right)+{\mathbf{Z}}_q\left(t\right)\end{array}$ (32)In these equations, $\mathit{\omega }\left(t\right)={\mathbf{J}}_r\left(\mathbf{q}\right)\stackrel{.}{\mathbf{q}}\left(t\right)$ ; $\mathrm{Log}:\text{SO}\left(3\right)\to {\mathbb{R}}^3$ denotes the matrix Logarithm Map (Lynch and Park, 2017; Murray et al., 1994); k _r , b _r denote the rotational stiffness and damping coefficients, respectively; R₀(t) ∈ SO(3) is the virtual orientation to which k _r and b _r are connected.

To generate a goal-directed discrete movement for both position and orientation of the robot’s end-effector, the motions for p₀(t) and R₀(t) to end at the goal position g and goal orientation R _g ∈ SO(3) are planned independently and simply added. Given a kinematically redundant robot, with positive k _r and b _r values and with a positive definite joint damping matrix B _q , p(t) → g and R(t) → R _g (Kim et al., 2011).

3.10.3. Simulation example

Consider a KUKA LBR iiwa robot, which has seven revolute joints (i.e., n = 7). We generated a single discrete movement for both position and orientation. The code script for the simulations can be found in main_redundant_pos_orientation.py

For position, a minimum-jerk trajectory was used (equation (17)). For orientation, a geodesic curve between initial R _i ∈ SO(3) and goal R _g orientations was used (Park and Ravani, 1994, 1997): ${\mathbf{R}}_{des}\left(t\right),{\mathbf{R}}_0\left(t\right)=\{\begin{array}{ll}{\mathbf{R}}_i\mathrm{Exp}\left(\frac{t}{T}\mathrm{Log}\left({\mathbf{R}}_i^{\mathrm{T}}{\mathbf{R}}_g\right)\right)& 0\le t<T\\ {\mathbf{R}}_g& T\le t\end{array}$ In this equation, $\mathrm{Exp}:{\mathbb{R}}^3\to \text{SO}\left(3\right)$ and $\mathrm{Log}:\text{SO}\left(3\right)\to {\mathbb{R}}^3$ are the matrix Exponential and Logarithm maps, respectively (Lynch and Park, 2017; Murray et al., 1994).

As shown in Figure 14, a goal-directed discrete movement for both position and orientation was achieved using both approaches. For DMPs, the approach used a feedback controller (equation (31)) with position and orientation trajectories generated by separate DMPs. For EDAs, the controller reused the impedance controller of Section 3.9.3 and simply superimposed an additional impedance operator Z _r for orientation. As demonstrated in Section 3.5 and Section 3.9, this example again emphasizes the simplicity due to the modular property of EDAs.

Figure 14.

Control of end-effector position and orientation while managing kinematic redundancy using Dynamic Movement Primitives (DMPs) and Elementary Dynamic Actions (EDAs). Successive frames of a MuJoCo model of a KUKA LBR iiwa14 robot using DMPs (A) and EDAs (B). The end-effector’s orientation R(t) is depicted with three orthogonal red, green, and blue arrows, while the desired orientation, R _des (t) (or R₀(t)) is depicted with three black arrows. Black dotted lines depict a minimum-jerk position trajectory (equation (17)), while the desired orientation followed a geodesic curve (Park and Ravani, 1994, 1997) from the start to the goal orientations. In detail, given the initial orientation R _i and goal orientation R _g , R _des (t) for DMPs (R₀(t) for EDAs) followed ${\mathbf{R}}_{des}\left(t\right)={\mathbf{R}}_i\mathrm{Exp}\left(t\mathrm{Log}\left({\mathbf{R}}_i^{\mathrm{T}}{\mathbf{R}}_g\right)\right)$ for t ∈ [0, 1], where Exp and Log are matrix Exponential and Logarithm maps, respectively (Lynch and Park, 2017; Murray et al., 1994). Compared to EDAs, control based on DMPs achieved perfect tracking for both position and orientation—the end-effector’s position and orientation perfectly overlapped the desired position and orientation, respectively. Despite the non-negligible tracking error for EDAs, the end-effector eventually converged to the goal position and orientation. As before, tracking error may be reduced by appropriate choice of impedance parameters. Parameters of DMPs for task-space position were identical to those in Figure 5. Parameters of DMPs for task-space orientation: α _z = 2000, β _z = 1000, N = 50, τ = 1.0, p = 100. Other controller parameters for DMP: Λ _p = 10I₃, Λ _r = 10I₃, Λ _q = 3I₇. Parameters of EDAs: K _p = 1600I₃ N/m, B _p = 800I₃ N⋅s/m, k _r = 80 N⋅m/rad, b _r = 8 N⋅m⋅s/rad, B _q = 20I₇ N⋅m⋅s/rad.

4. Combining the two approaches

In this section, we provide practical implementations that combine DMPs and EDAs on real robot hardware. The benefit of the combination is to exploit the modularity of EDAs and the robust trajectory generation capabilities of DMPs. Merging these two motor-primitive approaches not only emphasizes the individual strengths of each approach, but also creates an efficient framework for programming multi-task control, especially for the combined control of both task-space position and orientation (Section 4.1, 4.2).

Using a KUKA LBR iiwa14, two examples are provided: a drawing and erasing task (Section 4.1) and a modular Imitation Learning approach (Section 4.2). For both tasks, both position and orientation of the robot end-effector were controlled. To derive the torque command for the robot, the controller based on EDAs in Section 3.10.2 was used (equation (32)). To encode p₀(t), the DMP of Section 3.3.1 was used (equation (15)). For encoding R₀(t), the DMP of Section 3.10.1 was used (equation (29)).

For control, KUKA’s Fast Robot Interface (FRI) was employed. For both tasks, the built-in gravity compensation was activated. The Forward Kinematic Map to derive p, R, and the Jacobian matrices J _p (q), J _r (q) were calculated using the Exp[licit]^TM-FRI Library. ⁵ The C++ robot control code is available at https://github.com/mosesnah-shared/KUKA_app.

4.1. Drawing and erasing task

In this section, we provide an overview of the approach presented in Nah et al. (2023b), which considered a drawing and erasing task on a planar table. For further details on implementation, the readers are referred to Nah et al. (2023b).

The robot first drew a trajectory provided by demonstration. Once the demonstration was complete, the robot retraced the drawn trajectory with an additional oscillatory movement to thoroughly erase it.

This task can be executed by combining Imitation Learning of DMPs with the modularity of EDAs. Specifically, given a demonstrated trajectory p _des (t), Imitation Learning with a two-dimensional DMP is employed to learn this trajectory. Once this trajectory is learned, it is used for p₀(t) of EDAs (equation (32)).

After the trajectory was drawn, erasure was achieved by overlaying an oscillatory circular trajectory with radius r and angular velocity ω₀ onto the time-reversed path of p _des (t), described by the equation: $\begin{array}{l}{\mathbf{p}}_0\left(t\right)=\{\begin{array}{ll}{\mathbf{p}}_{des}\left(T-t\right)+{\mathbf{p}}_{0,osc}\left(t\right)& \hspace{1em}0\le t<T\\ {\mathbf{p}}_{0,osc}\left(t\right)& \hspace{1em}T\le t\end{array}\\ {\mathbf{p}}_{0,osc}\left(t\right)=r\left[\cos \left({\omega }_{0}t\right),\sin \left({\omega }_{0}t\right),0\right]\end{array}$

The result is shown in Figure 15, ⁶ which shows the whole process of the drawing and erasing tasks. By merging Imitation Learning with EDAs (Figure 15(A)–(C)), the drawing (Figure 15(D)) and erasing tasks (Figure 15(E)) were successfully achieved. The key to this approach is the combination of Imitation Learning with the modularity of EDA. The trajectory q _des (t) was separately learned using Imitation Learning and directly combined with an oscillation p_0,osc(t). Note in passing that, using EDAs, the robot remained stable despite physical contact with the planar table.

Figure 15.

The drawing and erasing task using a KUKA LBR iiwa. A green pen was used for drawing. (A) Data collection of human-demonstrated p _des (t) which we aim to draw. The end-effector trajectory along X- and Y-coordinates were collected. (B) Time t versus X-coordinate (black line) and Y-coordinate (purple line) of p _des (t) (top row), ${\stackrel{.}{\mathbf{p}}}_{des}\left(t\right)$ (middle row), ${\ddot{\mathbf{p}}}_{des}\left(t\right)$ (bottom row). With a sampling rate of 333 Hz, a first-order finite difference method was used to calculate velocity and acceleration (left column). These trajectories were Gaussian filtered (right column) using MATLAB’s smoothdata function with a time window of 165 ms. (C) The resulting trajectory p _des (t) (black dashed line) generated by Imitation learning. (D) The drawing task was achieved by setting p₀(t) as p _des (t). The black dashed line depicts p₀(t) (or p _des (t)), green line depicts p(t). (E) The erasing task was achieved by superimposing an oscillation p _osc (t) onto a time-reversed p _des (t). The green pen was replaced by a rectangular eraser. The green line depicts p(t). For (C, D, E), trajectories were plotted in MATLAB and overlapped onto the drawing/erasing table.

It is important to note that the Z-coordinate of p₀(t), p_0,z(t), was not determined through Imitation Learning. Instead, a suitable value for p_0,z(t) must be chosen to ensure contact with the table surface. If the table’s Z-coordinate is h _z and the pen attached to the robot’s end-effector points downwards (toward the −Z direction), p_0,z(t) should be set slightly below the table’s height, h _z , with a small positive offset ϵ, that is, p_0,z(t) = h _z − ϵ. This choice along with the inclusion of Z _r in τ _in (t) with a constant R₀ (equation (32)), allowed for the drawing process to be carried out while maintaining orientation.

4.2. Modular imitation learning

Combining the superposition principle of mechanical impedances (equation (10)) with Imitation learning offers a unique advantage in robot programming. It allows for the separate learning of position and orientation trajectories, that can then be added through the superposition of mechanical impedances. Hence, different motions learned in different spaces can be separately learned and seamlessly combined, which thereby enables modular Imitation Learning.

An example application for modular Imitation Learning is shown in Figure 16. ⁷ For task-space position, the robot learned to draw the letter “M,” using Imitation Learning with a three-dimensional DMP in task-space (equation (15)) (Figure 16(A)–(C)). For task-space orientation, the robot learned to “lift up and down” (as in lifting a cup to the lips and then putting it back down), using Imitation Learning with DMP employing unit quaternions (equation (29)) (Figure 16(D)–(F)). For both task-space position and orientation, Locally Weighted regression was used for Imitation Learning. By exploiting the superposition of mechanical impedances, the motion of drawing letter M while lifting up and down can be generated by a linear combination of both robot commands (Figure 16(G)). Note that while the qualitative behavior of both trajectories is preserved in combination, small tracking errors still remain.

Figure 16.

Modular Imitation Learning. Data for position and orientation were collected from real robot hardware (B, E) and visualized using Exp[licit]-MATLAB library (Lachner et al., 2023) (A, C, D, F, G). (A) A demonstrated trajectory p _des (t) for task-space position. p _des (t) followed an “M”-shaped trajectory. (B) Time t versus p _des (t), ${\stackrel{.}{\mathbf{p}}}_{des}\left(t\right)$ , ${\ddot{\mathbf{p}}}_{des}\left(t\right)$ with smoothing (equation (15)). (C) Successive frames of the robot movement using p₀(t) = p _des (t) with fixed R₀ (equation (32)). The robot was able to draw the letter M, which was provided by demonstration. (D) A demonstrated trajectory ${\overrightarrow{\mathbf{q}}}_{des}\left(t\right)$ for task-space orientation. ${\overrightarrow{\mathbf{q}}}_{des}\left(t\right)$ generated “a lifting up and down” motion. (E) Time t versus e _des (t), ${\stackrel{.}{\mathbf{e}}}_{des}\left(t\right)$ , ${\ddot{\mathbf{e}}}_{des}\left(t\right)$ with smoothing (equation (29)). (F) Successive frames of the robot movement using R₀(t) = R _des (t), where R _des (t) is converted from ${\overrightarrow{\mathbf{q}}}_{des}\left(t\right)$ . For this example, p₀ was fixed (equation (32)). (G) The combination of p₀(t) = p _des (t) and R₀(t) = R _des (t) (equation (32)). By exploiting the modularity of EDAs, combining the robot motion learned in different spaces—drawing letter M for task-space position, lifting up and down for task-space orientation—was achieved by simple linear addition of mechanical impedances. Parameters of DMPs for task-space position and orientation were identical to those in Figure 14. For both task-space position and orientation, Locally Weighted regression was used for Imitation Learning.

5. Discussion

In Section 3, we presented detailed implementations of both DMPs and EDAs to solve nine control tasks. Moreover in Section 4, we showed how these two methods can be combined to exploit the advantages of both approaches.

5.1. Similarities between the two approaches

DMPs and EDAs both stem from the idea of motor primitives. Hence, both approaches share the same principle—using motor primitives as fundamental building blocks to parameterize a controller. DMPs parameterize the controller with a canonical system, nonlinear forcing terms, and transformation systems (Section 2.1). EDAs parameterize the controller with submovements, oscillations, and mechanical impedances (Section 2.2).

Robot control based on motor primitives provides several advantages, and we presented nine control examples. First, by parameterizing the controller with motor primitives, the approaches are consistent with a high level of autonomy for generating dynamic robot behavior. Once triggered, the primitive behaviors “play out” without requiring intervention from higher levels of the control system.

As a result, the computational complexity of the control problem is reduced. For instance, we showed that DMPs can be scaled to multi-DOF systems by synchronizing a canonical system with multiple transformation systems (Section 2.1.4). With Locally Weighted (or linear least-squares) regression of Imitation Learning, learning new motor skills is reduced to calculating the best-fit weights of the nonlinear forcing terms. The best-fit weights are learned by simple matrix algebra, which is computationally efficient (Section 3.2, 3.3). In fact, it was reported that this computational efficiency of DMPs achieved control of a 30-DOF humanoid robot (Atkeson et al., 2000; Ijspeert et al., 2002, 2013; Schaal et al., 2007).

For EDAs, by parameterizing the controller with motor primitives, the process of acquiring and retaining complex motor skills is simplified by identifying a reduced set of parameters, for example, the initial and final positions of a submovement (Section 3.8, 3.9, 3.10). These computational advantages are particularly prominent in control tasks associated with high-DOF systems, for example, manipulation of flexible, high-dimensional objects (Nah et al., 2020, 2021; 2023a).

Motor primitives offer a modular framework for robot control. By treating motor primitives as basic modules, acquisition or generation of new motor skill occurs at the level of modules or their combination (d’Avella, 2016). Once the modules are learned, one can generate a new repertoire of movements by simply combining or reusing the learned modules. As shown above in several simulation examples, that is, obstacle avoidance (Section 3.5), combination of discrete and rhythmic movements (Section 3.7), sequencing discrete movements (Section 3.8), managing kinematic redundancy while controlling position (Section 3.9) and orientation (Section 3.10), these tasks were simplified by the modular properties of DMPs and EDAs. This modularity provides strong adaptability and flexibility for robot control, as learning new motor skills by combining or reusing learned modules is intuitively easier than learning “from scratch.” However, it is worth emphasizing that the details of the modular property are significantly different for the two approaches (Section 5.2.5).

5.2. Differences between the two approaches

5.2.1. The need for an inverse dynamics model

For torque-controlled robots, DMPs require an inverse dynamics model, whereas EDAs do not (Section 3.1). An inverse dynamics model can introduce practical challenges and constraints when applying DMPs to torque-controlled robots. In particular, acquiring accurate models of the robot can be challenging and time-consuming. Moreover, even if one acquires an accurate robot model, a low-gain joint PD controller must be additionally included to account for model uncertainties (Section 3.2, 3.3) and unexpected disturbances (Section 3.4).

However, the drawbacks associated with an inverse dynamics model can be dismissed for position control in joint-space. As a result, the application of DMPs in joint-space position-controlled robots is straightforward and efficient. On the other hand, EDAs are inappropriate for position-controlled robots, since “mechanical impedance” (Section 2.1.2) outputs (generalized) force, which is eventually commanded as a joint-torque input to the robot.

Nevertheless, joint-space position control presents its own challenges when compared to torque-controlled approaches. One obvious challenge is the problem of kinematic transformation between the robot’s generalized coordinates and task-space coordinates (Hogan, 2022). In fact, these challenges were highlighted in the examples presented, for example, the problem of inverse kinematics and kinematic singularity (Section 3.3) and managing kinematic redundancy (Section 3.9, 3.10). These problems necessitate additional methods, for example, Damped Least-squares inverse to manage kinematic singularity (Section 3.3) or feedback control methods such as sliding mode control to manage kinematic redundancy (Section 3.9, 3.10). Note that the eight control methods to manage kinematic redundancy presented by Nakanishi et al. (2008) (which include the latter approach) assume feedback control using torque-actuated robots, not position-controlled robots.

Perhaps more important, position control introduces the likelihood of instability for tasks involving contact and physical interaction. Robot control involving physical interaction requires controlling the interactive dynamics between the robot and environment (Hogan, 2022). For this, using position control turns out to be inadequate (De Santis et al., 2008). A position-actuated robot fails to provide the level of compliance needed to achieve safe physical interaction. Moreover, interactive dynamics cannot be directly regulated independent of the environment. Instead of position control, we considered control methods using (ideal) torque-actuated robots to manage contact and physical interaction (Section 3.4).

5.2.2. Tracking control

In the absence of uncertainties and external disturbances, in principle, DMPs can achieve perfect trajectory tracking, both in task-space and joint-space coordinates (Section 3.2, 3.3, 3.6, 3.7, 3.9, 3.10). Using Imitation Learning (Section 2.1.4), tracking a trajectory of arbitrary complexity can be achieved. DMPs also allow online trajectory modulation of the learned trajectory, which was shown in the obstacle avoidance example (Section 3.5).

EDAs control neither position nor force directly (Section 2.2.4). Hence, non-negligible tracking errors arise unless high values of mechanical impedance are employed (Section 3.3, 3.6). One can reduce the tracking error by using higher stiffness and/or damping values (Figures 4 and 6). Nevertheless, a non-negligible tracking error still exists. To achieve higher accuracy for tracking control with EDAs for a given desired trajectory, an additional method to derive the corresponding virtual trajectory should be employed. For instance, exploring trajectory optimization methods to calculate the time course of impedances and virtual trajectories that produce the desired trajectory is an opportunity for future research.

In fact, the topic of merging trajectory optimization with EDA is only one example of a potentially rich research field, which is to integrate the mathematical framework of optimal control theory into EDA. The key point of EDA is to choose appropriate values of mechanical impedance and the virtual trajectory to achieve the control task. The parameters of EDA used for the nine control examples, especially for the mechanical impedances, were mostly chosen based on trial-and-error, and these parameters were sufficient to achieve the task. Nevertheless, other work advances EDA by optimizing the parameters based on specific cost functions, for example, tracking error for position (Buchli et al., 2011).

The presented examples mainly focused on tracking control of point-to-point goal-directed discrete movements. For EDA, the movements can be generated by submovements and their linear composition. For DMP, the movements can be achieved with or without the nonlinear forcing term, where the latter follows the response of a stable second-order linear system. However, if one considers path-constrained tracking control, the benefit of DMP’s Imitation Learning is more emphasized. Compared to EDA, where the path should be decomposed into multiple submovements and their combination (Rohrer et al., 2002, 2004), DMP can directly learn the path by Imitation Learning which can be used for path-constrained tracking control. The path learned by Imitation Learning can also be generalized by spatial and temporal scaling, which makes DMP favorable over conventional spline methods (Ijspeert et al., 2013). In fact, this specific benefit of DMP for path generation was combined with EDA to get the best of both motor-primitives approaches, as shown in the two examples presented in Section 4. In summary, while both EDA and DMP offer methods for generating goal-directed movements, DMP’s Imitation Learning presents a significant advantage in path-constrained tracking control, providing a more efficient and adaptable solution compared to EDA.

5.2.3. Contact and physical interaction

To guarantee robustness against contact and physical interaction, DMPs superimpose a joint-space PD controller on the feedforward torque command from the inverse dynamics model (Section 3.1). On the other hand, EDAs include mechanical impedances as a distinct class of primitives (Section 2.2.3). With an appropriate choice of mechanical impedance, the approach is robust against uncertainty and unexpected physical contact (Hogan, 2022). The dynamics of physical interaction can be directly controlled by modulating mechanical impedance (Section 3.4). By superimposing passive mechanical impedances, passivity is preserved (Section 3.9, 3.10, 4.1).

Note that the equation of PD control used for DMPs (equation (18)) is identical to a first-order joint-space impedance controller (equation (13)). However, care is required: they are identical only if the robot actuators are ideal torque sources and the impedance is specified in joint space. Impedance control is more general than PD control and not limited to first-order joint-space behavior; by definition, mechanical impedance determines the dynamics of physical interaction at an interaction port, which may, in principle, be at any point(s) on the robot (Won et al., 1997).

5.2.4. Managing kinematic singularity and redundancy

For DMPs, control in task-space requires solving an inverse kinematics problem (Section 3.3). This introduces the challenges of managing kinematic singularity and kinematic redundancy. The latter can be resolved by using any of many feedback control methods presented by (Ceccarelli, 2008; Nakanishi et al., 2008). Nevertheless, this requires feedback control based on an error signal. This introduces a non-negligible error, and instead of perfect tracking, asymptotic convergence is achieved. Moreover, these methods still involve a Jacobian (pseudo-)inverse. Consequently, an additional method to handle kinematic singularity should be employed. Importantly, null-space projection methods violate passivity (Section 3.9.1, 3.10.1), and advanced methods to guarantee the robot’s stability might be needed (Dietrich et al., 2015; Lachner, 2022).

For EDAs, explicitly solving the inverse kinematics is not required (Hogan, 1987). Seamless operation into and out of kinematic singularities is possible (Section 3.3.3). EDAs superimpose multiple mechanical impedances to manage kinematic redundancy (Section 3.9.3, 3.10.3). Unlike null-space projection methods, passivity is preserved (Hogan, 2022; Lachner, 2022).

5.2.5. Modularity in robot control

As discussed in the examples of Section 3 and Section 5.1, for both approaches, motor primitives provide a modular control framework for robot control, which thereby simplifies programming of multiple control tasks (Section 3.5, 3.7, 3.8, 3.9, 3.10). Nevertheless, it is important to note that the extent of modularity and its practical implications significantly differ between these two approaches.

A clear distinction is evident when combining multiple movements. For DMPs, by design, discrete and rhythmic movements are generated by different DMPs. Hence, discrete and rhythmic movements cannot be simply superimposed (Section 3.7). Moreover, to sequence discrete movements, the goal location g(t) of the previous discrete movement is modulated (equation (26)) (Section 3.8).

On the other hand, for EDAs, sequencing and/or combining multiple movements can be seamlessly conducted at the level of the virtual trajectory simply by adding components (Section 2.2.4). Recall that a combination of discrete and rhythmic movements was achieved by simply adding submovements and oscillations (Section 3.7). For sequencing discrete movements, the subsequent discrete movement was superimposed without modifying the previous movement (Section 3.8). These properties provide a notable degree of simplicity and modularity, as individual motions can be separately planned and simply added without further modification.

The superposition principle of mechanical impedances enables breaking down complex tasks into simpler sub-tasks, solving each sub-task with a specific module, and simply adding the outputs of these modules to solve the original problem (Section 3.5). Using mechanical impedances with this “divide-and-conquer” strategy, the overall complexity of the control problem can be significantly reduced. A task may be achieved by simply reusing the impedance controllers of component tasks without modification. This modular property of EDAs is in contrast with DMPs. While the learned weights of the nonlinear forcing terms can be reused for a single DMP, multiple DMPs cannot be simply combined by merging the learned weights of different DMPs.

5.3. Combining the best of both approaches

Despite these differences, both approaches can be combined to leverage their respective advantages and alleviate their limitations. EDAs are robust against uncertainty and have advantages for physical interaction. DMPs can easily learn and track trajectories of arbitrary complexity. A combination of both approaches allows robustness against uncertainty and physical interaction, while enabling efficient learning and tracking of trajectories with arbitrary complexity (Table 2).

Table 2.

The Advantages of DMP and EDA That Are Combined to Get the Best of Both Approaches.

	Dynamic movement primitives (DMPs)	Elementary dynamic actions (EDAs)
Advantages	Capable of generating a wide range of trajectories using Imitation Learning	Modularity at the level of motion planning, interactive dynamics, and robot command
Combining the best of both approaches	Generating the virtual trajectories of EDA using Imitation Learning of DMP. Section 4.1: Merging discrete DMP with oscillation for the drawing and erasing task. Section 4.2: Separately learn trajectories for task-space position and orientation, and then linearly combine at the level of robot command

In Section 4, we provided two examples that demonstrated the benefits of combining both approaches. The key idea was to use Imitation Learning of DMPs to generate the virtual trajectories of EDAs. With this, the drawing and erasing task was successfully achieved (Section 4.1). Thanks to EDAs favorable stability property, the robotic manipulator remained stable despite contact and physical interaction. With Imitation Learning of DMPs and the kinematic modularity of EDAs, the erasing motion was achieved by simply superimposing an oscillation onto the (reversed) drawing motion. Moreover, modular Imitation Learning was achieved where different motions learned in different spaces were separately learned and linearly combined (Section 4.2).

Along with Section 4, further examples illustrate combinations of the two approaches. Superimposing a low-gain PD controller onto the feedforward torque command of DMPs (Section 3.4) can be regarded as an example of combining both approaches. Given an ideal torque-actuated robot, a first-order joint-space impedance controller (an example of an EDA) is added to a feedforward torque command based on DMPs. Robustness against uncertainty or physical contact can also be achieved by superimposing other mechanical impedances, for example, a first-order position task-space impedance controller (equation (16)). Additionally, for tasks combining discrete and rhythmic movements (Section 3.7) and sequencing discrete movements (Section 3.8), one can instead encode g(t) using submovements and/or oscillations of EDAs rather than a first- (equation (25)) or second-order linear systems (equation (26)). This further provides flexibility of motion planning for DMPs, and the motion need not be confined to the response of a low-order linear system.

6. Conclusion

In this paper, we provided a detailed comparison of two motor-primitive approaches in robotics: DMPs and EDAs. Both approaches utilize motor primitives as fundamental building blocks to parameterize a controller, enabling highly dynamic robot behavior with minimal high-level intervention.

Despite this similarity, there are notable differences in their implementation. Using simulation, we delineated the differences between DMPs and EDAs through nine robot control examples. While DMPs can easily learn and track trajectories of arbitrary complexity, EDAs are robust against uncertainty and have advantages for physical interaction. Accounting for the similarities and differences of both approaches, we provided real robot implementation of how DMPs and EDAs can be combined to achieve a rich repertoire of movements that is also robust against uncertainty and physical interaction.

In conclusion, control approaches based on DMPs, EDAs or their combination offer valuable techniques to generate dynamic robot behavior. By understanding their similarities and differences, researchers may make informed decisions to select the most suitable approach for specific robot tasks and applications.

Supplemental Material

Supplemental Material