Invariant Descriptors of Motion and Force Trajectories for Interpreting Object Manipulation Tasks in Contact

Invariant descriptors of point and rigid-body motion trajectories have been proposed in the past as representative task models for motion recognition and generalization. Currently, no invariant descriptor exists for representing force trajectories, which appear in contact tasks. This article introduces invariant descriptors for force trajectories by exploiting the duality between motion and force. Two types of invariant descriptors are presented depending on whether the trajectories consist of screw or vector coordinates. Methods and software are provided for robustly calculating the invariant descriptors from noisy measurements using optimal control. Using experimental human demonstrations of 3-D contour following and peg-on-hole alignment tasks, invariant descriptors are shown to result in task representations that do not depend on the calibration of reference frames or sensor locations. The tuning process for the optimal control problems is shown to be fast and intuitive. Similar to motions in free space, the proposed invariant descriptors for motion and force trajectories may prove useful for the recognition and generalization of constrained motions, such as during object manipulation in contact.


I. INTRODUCTION
I N human-robot interaction (HRI) there is a need for generalizable task models that allow robots to reactively adapt the robot's execution and to robustly recognize human actions in continuously changing environments.For execution purposes, these models should be easily adaptable to deal with obstacles, different starting/ending/via points, robot platform constraints, and human-robot interaction.For recognition purposes, these models should be sufficiently general to recognize the human's action at different locations in space, from different viewpoints, and with different execution styles (e.g.varying scale and motion profile).The focus of this paper is on extracting task representations that involve the manipulation of objects in contact, where not only the relative motion between objects but also interaction forces are relevant as input for the task modeling process.Typical applications include surface following tasks and assembly tasks.
Task models in HRI are commonly learned from one or multiple human demonstrations of the task using Learning by Demonstration [1].The challenge is to achieve a representative task model without the need for a high number of demonstrations.
One approach to achieve representative task models from few demonstrations is to exploit invariance in the task.Here, invariance refers to properties of the task that remain unchanged under certain transformations such as changes in reference frame.Invariance has proven useful in various fields of robotics.Typically, the objectives were: (1) to ensure that exactly the same task execution or behavior of a robot system was obtained within its workspace, and (2) to ensure that the quantitative characterization of such robot action using a metric resulted in exactly the same number(s) regardless of where and how the robot action was recorded or measured.Below, some examples of the importance and use of invariant solutions in different robotics contexts are briefly reported.
For the control of free-space robot motion such as mobile robots and UAVs, invariance was achieved by representing the trajectory tracking error in such a way that the same closed-loop behavior for trajectory tracking and disturbance rejection was obtained regardless of the direction in which the robot was moving [2].For constrained motion, i.e., motion in contact, which is the subject of this paper, the risk of obtaining non-invariant control solutions is more pronounced since both motion (translation and rotation) and interaction wrenches (i.e.forces and moments) have to be controlled simultaneously.For example, as pointed out in [3], it does not make sense to define orthogonality between two twists (i.e.translational and rotational velocities), or between two contact wrenches, because of (1) dimensional inconsistency, (2) dependence on the choice of units, and (3) dependence on the choice of the origin of the coordinates.In contrast, reciprocity between a twist and a wrench is an invariant concept.In [4], [5] a firm geometric foundation for kinestatics and invariant hybrid force/position control was laid based on screw theory.Taking this into consideration, invariant solutions were proposed, for example to obtain dynamic models of redundant and constrained robot manipulators [6], or in applications of hybrid force/position control [7].
In estimation, invariance refers to obtaining the same values for the estimated quantities regardless of the reference frame in which the estimator's states are represented.For example, the Invariant Extended Kalman Filter [8] formulated the state correction term in an invariant way using Lie group theory and showed a better convergence of the estimation compared to a standard Extended Kalman filter with a linear correction term.The Invariant Extended Kalman Filter was applied, for instance, to attitude estimation for UAVs [9].
In motion generation, invariance refers to obtaining exactly the same physical desired motion trajectory regardless of the reference frame in which the trajectory coordinates are expressed.Human reaching motions, represented as point trajectories, were found to correspond well to a minimum-jerk profile, an invariant property that was shown to be equivalent to having a constant equi-affine curvature [10].Based on this property, new trajectories can be generated or adapted from reference trajectories in an invariant way [11].This and other invariant properties such as shape preservation were considered for the deformation of trajectories using global transformations [12].For rigid-body motion, approaches were proposed where the trajectories are either generated in a leftinvariant way, meaning independent of the world reference frame, or in a right-invariant way, meaning independent of the body reference frame [13].This was also taken into account for the adaptation of reference or demonstrated rigid-body trajectories [14], [15].
Invariance in action recognition refers to recognizing actions regardless of where and how they are recorded.This can be done by comparing the recorded action instance and its model using invariant distance metrics.For example, local invariant signatures were proposed for describing and recognizing object outlines [16] and trajectories [17].These invariant representations were shown to be invariant to occlusion, rotation, translation, and scaling.
These previous works showed that invariant representations offer many advantages that are relevant for different robotics applications.In the remainder of the paper we focus on invariant descriptors for trajectories.These descriptors can be used to build representative invariant models of a demonstrated task in a human-robot interaction context with fewer demonstrations compared to conventional methods using coordinatebased descriptors.In addition, thanks to their invariant properties, these models enable robust action recognition in different recording contexts and exhibit good generalization properties in motion generation.

A. Invariant trajectory descriptors for motion
Invariant trajectory descriptors are features of the trajectory that remain unchanged under certain transformations of the original trajectory coordinates, such as changes in reference frame, scale, and time profile.Invariant descriptors can be categorized in many types [18] depending on the features that are extracted from the curve.Many have been proposed, both for point trajectories [16], [19]- [22] and for rigid-body trajectories [23]- [27].
This paper focuses on a particular type of invariant descriptor, referred to as differential invariants [16], [18].These invariants correspond to differential-geometric properties of the trajectory and provide a local invariant representation that is robust to occlusion and segmentation errors.Furthermore, they allow to reconstruct and generate new trajectories from the invariant descriptor by integrating corresponding differential equations.The definition of such local differential-geometric invariants can be seen as an application of Cartan's method of moving frames [28], which consists of defining a natural moving frame in the considered manifold, after which the derivative of the frame is shown to correspond to a set of geometric invariants.
For point motion in Euclidean space, the differential invariants correspond to curvature and torsion defined in a moving frame referred to as the Frenet-Serret frame.Originally introduced in vision for characterizing and recognizing objects by their outlines [16], [29], the same concepts of curvature and torsion were applied for describing and recognizing point motion trajectories [17], [20].
For rigid-body motion, differential invariants were introduced based on the concept of the screw axis in screw theory [25].A moving frame was defined on the screw axis of which the motion was described by a set of invariants referred to as screw invariants.These invariants are bi-invariant, meaning independent of both the world reference frame and the reference frame attached to the body.An alternative invariant descriptor was proposed by defining two separate Frenet-Serret frames for the position and orientation of the rigid body [26], [27], resulting in a simpler set of formulas.However, invariance for the choice of reference point for translation (the body frame's origin) was lost.
Invariant trajectory descriptors for motion have, for the large majority, been applied in motion recognition [20], [26], [27].Trajectory generation using invariant descriptors has also been explored.The invariant descriptors are capable of both reconstructing the original trajectory [25], as well as reproducing the trajectory at new locations and in different directions by applying the correct initial values when integrating the differential equations [30].Building further on that, trajectory adaptation methods based on optimal control and invariant descriptors were proposed [15].The goal was to adapt a reference trajectory to comply with new constraints while minimizing deviation from the invariant descriptor of the reference trajectory.This was shown to result in excellent extrapolation capabilities [15].
Calculating differential invariants from measurement data is challenging due to noise-sensitivity and singularities.Singularities are defined as instances along the trajectory where components of the descriptor are arbitrarily defined.Near these singularities, small variations in the trajectory (e.g., due to measurement noise or irrelevant human variations) will result in large variations in the ill-defined components of the descriptor.To remedy this, efforts have been spent to develop robust approaches to calculate invariant descriptors, ranging from discrete approximations of the descriptor [16], [17], [20], [26] to averaging the descriptor over segments of the curve [31].More recently, approaches based on optimal control were introduced for calculating point trajectory descriptors [32] and rigid-body trajectory descriptors [33].Using optimal control, the local descriptor is calculated in a window of measurements while regularization terms deal with singularities and noise by smoothing the invariants in a local neighborhood.A

Trajectory
Applied to motion Applied to force vectors c disadvantage of current optimal control-based approaches is their dependency on many tuning parameters for weighting the trajectory accuracy cost and regularization costs, and their long calculation times.

B. Paper objective and contributions
Up till now, invariant trajectory descriptors for interaction forces with similar benefits as the invariant descriptors for motion have not been proposed.This paper's main objective is to extend the existing concepts and methodology for motion invariants towards interaction forces by exploiting duality between motion and force.The resulting invariant descriptors can be used to model motion and force trajectories that occur in motion-in-contact tasks.Table I summarizes this duality for the two cases where motion and force trajectories are either represented with vector coordinates or with screw coordinates.
The main contribution of this paper is the introduction of two types of local differential invariant descriptors for force and moment trajectories.The first type is referred to as the vector invariant descriptor for vector trajectories and is applied to force and moment vectors separately.The second type is referred to as the screw axis invariant descriptor for screw trajectories and is applied to a screw wrench, combining force and moment in one entity.The duality is shown between the new invariant descriptors for force and moment and the existing invariant descriptors for motion: the screw invariants for motion [25] and the vector invariants for translation and rotation [26], [27].Furthermore, we explain how to calculate the invariant descriptors for motion and force using numerical optimal control schemes, which filter out the effects of measurement noise and improve the behavior of the descriptor near singularities.Triggered by the higher noise levels of measured force and moment trajectories compared to measured motion trajectories, we improved the formulation of these numerical optimal control schemes compared to [33], resulting in much more intuitive and faster tuning as well as increased robustness of the calculation.The invariant properties of the proposed descriptors were experimentally verified using a human-demonstrated 3D contour following task and a peg-on-hole alignment task.All software used in the validation experiment, including data, has been made publicly available [34].
The remainder of the paper is structured as follows.Section II provides the necessary background on rigid-body kinematics and statics.Sections III and IV introduce local invariant 6×1 screw twist of rigid body with the first three elements the rotational velocity vector aω and the last three elements the translational velocity vector av a aw a 6×1 screw wrench acting on the rigid body with the first three elements the force vector af and the last three elements the moment vector am a b a S 6×6 screw transformation matrix that, when multiplied with a screw twist or wrench, transforms the screw from {b} to {a} descriptors in a general form, resulting in existing descriptors when applied to motion trajectories and new descriptors when applied to force and moment trajectories.Sections III and IV deal with vector and screw invariant descriptors, respectively.Section V highlights a practical 3D contour following task, while Section VI examines a peg-on-hole alignment task.
Section VII provides a discussion about the benefits and limitations of the proposed concepts and methodology, it points at future work and presents a conclusion.

II. PRELIMINARIES
This section reviews essential background and introduces notation as summarized in Table II.

A. Rigid-body kinematics
The motion trajectory of a rigid body is commonly represented by attaching a reference frame {b} to the body and expressing the position and orientation of this frame as a function of time t with respect to a fixed reference frame {w}, also referred to as the world reference frame.For all variables introduced below, the explicit dependency on time t is omitted to simplify the expressions.
The position is represented by the 3-dimensional displacement vector b w p from the origin of {w} to the origin of {b}.The orientation is represented by the rotation matrix b w R, which expresses the coordinates of the unit vectors of the moving body frame {b} with respect to the world frame {w}: b w R = e x e y e z . ( Position and orientation can be combined into the homogeneous transformation matrix or pose matrix b w T : the fixed reference frame {w} in the world and the reference frame {b} on the rigid body.
The velocity of a rigid body can be fully characterized by two 3-dimensional vectors: the rotational velocity vector ω and the translational velocity vector v of the origin of the considered reference frame.These vectors are summarized with a 6-dimensional screw [35], the twist t = ω v .
To define the time-derivative of the homogeneous transformation matrix, Ṫ = dT dt , we first introduce the operator [•] × .When applied to the vector ω = [ω x ω y ω z ] T , it results in a 3 × 3 skew-symmetric matrix: and when applied to a screw t, it results in the following 4×4 twist matrix [t] × , which is an element of the Lie algebra se(3): The definition of the time-derivative of the homogeneous transformation matrix Ṫ depends on the choice of reference frame in which the twist t is expressed.Choosing the body frame {b}, the body twist b t b contains the rotational velocity vector b ω expressed in {b} and the translational velocity vector b v b of the origin of {b} and expressed in {b} (see Fig. 1a).The derivative of the pose matrix b w Ṫ is then defined as: The twist b t b has the property of being invariant for a change of {w} (left-invariance).Alternatively, choosing the world frame {w}, the spatial twist w t w consists of the rotational velocity vector w ω expressed in {w}, and the translational velocity vector w v w of the point on the rigid body that instantaneously coincides with the origin of {w} expressed in {w}.The pose derivative b w Ṫ is then defined as: where [ w t w ] × is now multiplied on the left of b w T since the twist is expressed in {w}.The twist w t w has the property of being invariant for a change of {b} (right-invariance).
The relation between twists w t w and b t b is provided by the 6 × 6 screw transformation matrix S( b w T ) (also referred to as the adjoint matrix [36]), which is constructed from elements of the pose matrix b w T : Twists b t b and w t w represent first-order kinematics of the trajectory and are invariant to the choice of world frame and body frame, respectively.But neither w t w nor b t b is invariant for the choice of both reference frames {w} and {b}.Trajectory descriptors that are invariant to both reference frames (also referred to as coordinate-invariance) will be discussed in Sections III and IV.

B. Duality between twist and wrench as general screws
All forces and moments acting on a rigid body can be reduced to a single resultant force vector and a single resultant moment vector.One option is to express these vectors in the body frame: b f represents the resultant force expressed in {b}, and b m b represents the moment with respect to the origin of {b} and with coordinates expressed in {b} (see Fig. 1b).Similar to the twist, a compact notation for the resulting force and moment is the 6-dimensional screw wrench: Transforming the wrench from the body frame {b} to the world frame {w} is done using the same 6 × 6 screw transformation matrix S( b w T ) as for the twist: This analogy between twists and wrenches is also referred to as the duality between velocity and force.Both twist and wrench can be interpreted as a 6-dimensional screw [35], meaning they can always be represented by a vector lying along a line in space and the moment of that vector around the line.

C. Time-invariance using a geometric progress variable
In the previous subsections, motion and force trajectories were assumed to be a function of time t.To obtain a timeinvariant trajectory representation, i.e., a trajectory representation that is invariant with respect to the motion profile, a geometric progress variable ξ must be chosen.
For tasks where translation is dominant, such as in contour following, the arc length is a natural choice.The arc length is calculated by integrating the norm of the translational velocity vector v over time.For tasks where rotation is dominant, such as in opening a door, a natural progress variable would be the integral of the norm of the rotational velocity vector ω.If translation and rotation are both important, then a progress parameter can be defined combining the translation and rotation, for example, using a linear combination [25]: where w weights the relative contribution.In [25], further details are provided for obtaining a progress value that is invariant for the choice of reference point on the body for the translation v.
For the pose T and wrench w, the reparameterization from time t to progress ξ is straightforward: T (ξ) = T (t(ξ)) and w(ξ) = w(t(ξ)).(11) In practice, since measurement data is often only available in a discretely sampled form, the reparameterization is achieved through numerical interpolation.
For the velocity, the reparameterization from time t to progress ξ requires an additional step since velocities represent derivatives of the trajectory.As an example, the reparameterization of the translational velocity vector v, as the derivative of the position p, can be worked out as: where the chain rule of differentiation was applied in the second equality to obtain the relation with time t, and the progress rate ξ = dξ dt is the derivative of the progress with respect to time.Similar derivations can be made for the rotational velocity vector ω and the twist t: Equations ( 11)-( 13) can be used in the other direction to apply a motion profile ξ(t) and obtain the pose, wrench, translational velocity, rotational velocity or twist as a function of time.This is relevant, for example, when a trajectory was planned in the progress domain ξ and needs to be executed with a specific timing ξ(t).A special choice for the progress variable is one where ξ = 0 at the beginning of the trajectory and ξ = 1 at the end of the trajectory.This results in a reparameterization that is independent of the scale of the trajectory [25].

III. INVARIANT DESCRIPTOR FOR VECTOR TRAJECTORIES
This section derives a local invariant descriptor for a general vector trajectory c(ξ).The trajectory c(ξ) is assumed to be time-invariant, i.e. parameterized as a function of a geometric progress ξ as discussed in Section II-C.To keep the formulas concise, the explicit dependency on ξ is omitted in the following formulas.The derivative of c with respect to the progress ξ is represented using the prime notation: c ′ = dc dξ .

A. Definition of invariants
The vector invariants are based on a generalization of the Frenet-Serret equations for translation [37], where c corresponded to the translational velocity vector.While in [26], [27] the Frenet-Serret equations were extended to rotational velocity trajectories, the aim here is to exploit the duality between velocity and force so that the invariants can be presented in a general form, applicable to rotational and translational velocity trajectories, as well as to force and moment trajectories.
The vector invariants are defined and expressed in a local moving frame, referred to as the Frenet-Serret frame (FS frame), which is constructed as follows.Its first axis, the tangent e t , corresponds to the unit vector in the direction of c.The second axis corresponds to the normal e n of the trajectory.The third axis, the binormal e b , follows directly from the cross-product of e t and e n .These unit vectors are defined as follows from c and c ′1 : In differential geometry, the linear span of e t and e n is referred to as the osculating plane, i.e. the plane in which the vector c rotates instantaneously, while e b represents the axis about which this rotation occurs.The complete orientation of the FS frame is expressed using a rotation matrix R: R = e t e n e b , where the tilde ∼ indicates it is a local moving frame, in this case the FS frame.See Fig. 2 for a visualization of the frame.
The invariants are then defined as follows.The first invariant c is the magnitude of the vector c along the tangent: This invariant is referred to as the object invariant since it relates to the object.For example, c can represent the magnitude of the translational or rotational velocity vector of the object, or the magnitude of the force or moment vector applied to the object.The vector c can be retrieved from the object invariant c using: The second and third invariants ω κ and ω τ govern the firstorder kinematics of the local FS frame2 : where ω κ is known as the curvature rate and ω τ as the torsion rate, with the following explicit formulas [37]: For a naturally parameterized point trajectory, with c the point's velocity, ω κ and ω τ correspond to the local curvature κ and torsion τ of the point's trajectory.Equation ( 18) can be written compactly using (3): where i is referred to as the moving frame invariant since it models the first-order kinematics of the moving frame.
Only the orientation is relevant for the FS frame, not its position.Therefore it is sometimes referred to as an orientation frame [38].For visualization purposes, its evolution is often shown by representing it at the point on the motion trajectory associated with the progress variable ξ.This was done in Fig. 2, which shows the definition of the moving frame and the corresponding invariants for a point trajectory.
Singularities: Singularities are defined as instances where axes of the moving frame and the corresponding invariants are not uniquely defined anymore.A first type of singularity is when c = 0.The object invariant c is then zero, while the moving frame invariants ω κ and ω τ are undefined.In other words, the complete moving frame is arbitrary.
A second type of singularity is when the vector c has an unchanging orientation at some instance (c × c ′ = 0), remaining parallel to itself.The first moving frame invariant ω κ is then zero while the second moving frame invariant ω τ is undefined.In other words, the tangent of the moving frame is well-defined, but the normal and binormal axes are arbitrary.
Invariant properties: If the trajectory c(ξ) is expressed as a function of a geometric progress variable ξ as explained in Section II-C, then the descriptor (c, i) is also time-invariant.Furthermore, the descriptor is invariant with respect to the choice of world frame {w} (both origin and orientation) in which the coordinates of the given vector trajectory c(ξ) are expressed.It is evident that a change of {w} merely results in the same change of orientation for the vectors c, c ′ , and c ′′ , resulting in the same object invariant c, the same FS frame for all ξ, and hence also the same moving frame invariant i.If the coordinates of c are instead expressed in {b}, the descriptor is still invariant for changes in orientation of the body frame {b}.However, not all considered instances of the descriptor are invariant with respect to the origin of body frame {b}.As pointed out in [27], the translational velocity vector v depends on the reference point on the moving object chosen to express its velocity, hence on the origin of {b}.Similarly, the moment vector m depends on the reference point for expressing the moment acting on the object, hence on the origin of {b}.The dependency of v and m, and hence also of their invariant descriptors, on the origin of {b} is a disadvantage compared to the screw invariants discussed in Section IV, which do not exhibit this dependency on the chosen reference point on the moving body.
Trajectory reconstruction: Reconstructing the trajectory c(ξ) from the invariant descriptor ( c(ξ), i(ξ)) requires integrating the differential equation ( 20) from an initial value of the moving frame R. In general, there exists no closed-form solution for the integration.Instead, the problem is discretized so that R k+1 at sample k + 1 is found from R k at sample k as follows: where we assume that i k remains constant over the integration step ∆ξ.The matrix exponential operator exp(•) maps the skew-symmetric matrix [i k ] × ∆ξ, which is part of the Lie algebra so(3), into the corresponding change in orientation, which is part of the Lie group SO( 3).An explicit expression of the matrix exponential can be formulated using Rodrigues' rotation formula [36].Given the initial value R 1 at k = 1, the moving frame can be reconstructed for all remaining samples k = 2...N .From the moving frames R k , the corresponding vector trajectory c k at k : 1...N is found using (17): Equations ( 21) and ( 22) are referred to as the moving frame equation and object trajectory equation, respectively.Table III on page 10 provides an overview of the vector invariants framework with application to translational velocity v, rotational velocity ω, force f , and moment m trajectories.

B. Robust calculation of invariants using optimal control
Using the analytical formulas to calculate the vector invariants c and i has three pitfalls: (1) near singularities, some of the invariants and axes of the moving frame R are illdefined, hence the behavior of c and i is not well-captured in such case; (2) the high-order derivatives cause sensitivity to noise, especially near singularities; and (3) the open-loop reconstructed trajectory using ( 21)-( 22) is not guaranteed to match the original trajectory due to integration drift.
In addition, the invariants c and ω κ are always positive by definition from ( 16)- (19) resulting in frame flipping issues: the FS frame must flip 180 • when the direction of c reverses or when the curve traced by the vector c has an inflection point.Other definitions of moving frames exist to minimize the frame's rotation, such as Bishop frames [39], but these frames are not locally defined anymore.
To deal with these problems, we formulate an Optimal Control Problem (OCP) to calculate the vector invariants c and i in a robust way.An OCP is a type of optimization problem that, for a given dynamical system, tries to find control inputs that minimize a specific cost function over a horizon.Here, the objective of the optimization problem is to find vector invariants that reconstruct the given measured vector trajectory.The latter can be quantified by introducing a trajectory reconstruction cost ∆c MSE containing the meansquared error (MSE) between the reconstructed trajectory c k and measured trajectory c meas k for k = 1...N : To deal with measurement noise and the effects of singularities, a regularization cost is introduced.In general, the regularization cost can be any function of invariants, including derivatives of invariants.In [33], the terms in the regularization cost were chosen as the derivatives of all invariants and the absolute values of the moving frame invariants.The trajectory reconstruction cost and regularization cost were added in the objective function of the OCP, and as a result they had to be weighted with appropriate weighting factors, requiring an extensive tuning process.
To reduce tuning efforts and improve interpretability of chosen parameters, this paper proposes a new formulation of the OCP.The trajectory reconstruction cost is moved to the constraints of the OCP and is bounded by a chosen tolerance value ϵ.For vector trajectories, this is worked out as follows: in which ϵ c is interpreted as the desired tolerance on the trajectory reconstruction cost over the entire horizon.The subscript in ϵ c refers to the type of trajectory that is considered.The value of ϵ can be set by the user based on the expected accuracy of the measurement system or the tolerable reconstruction error for the given application.
For the regularization term in the objective function, we take the squared norm of the moving frame invariants i.Since i signifies the rotational change of the moving frame, this will result in a more stable and smooth evolution of the moving frame R. To avoid frame flipping issues, all invariants are allowed to become positive or negative, resulting in a smoother evolution of the FS frame.
The complete OCP to calculate the vector invariants c and i from a given sequence of N measured vectors c meas k is: minimize subject to: Besides the trajectory error constraint (24), the constraints consist of the trajectory reconstruction equations ( 21)-( 22), while constraint (26) imposes orthonormality conditions on the initial moving frame.The orthonormality of R is preserved throughout the entire horizon through the matrix exponential integrator in (21).
Initialization: The OCP ( 25)-( 26) is a non-convex optimization problem and therefore multiple local minima may exist.The solution will converge towards one of the local minima depending on the initialization of the unknown variables.To obtain invariant solutions regardless of the reference frame in which the variables are represented, the variables must be initialized in an invariant way.One possibility is to calculate initial values for the invariants and moving frames using the analytical formulas ( 15)-( 19), as was done in [33].However, this initialization then suffers from the same three pitfalls mentioned at the beginning of this subsection.
We propose an alternative initialization approach here based on calculating an average moving frame, first introduced in [40] for screw trajectories, but adapted here to vector trajectories in three steps.First, calculate the covariance matrix C c associated with the measured vector trajectory c meas k : Second, the orthogonal matrix of eigenvectors U c is calculated from C c .This matrix U c is invariant for changes in reference frame by definition and can be interpreted as a sort of average FS frame for the whole trajectory.The signs of the first two eigenvectors are chosen unambiguously by making the average value of the measurement vector c meas k positive along these directions.The third axis then follows from the cross product of the first two axes.All moving frames R k are initialized with the average moving frame U c .
Third, to be consistent with the average moving frame, the moving frame invariants ω κ and ω τ are initialized with zero values.The object invariant c k is initialized by projecting the vector trajectory c meas k onto the average moving frame and selecting the component along the first axis for all samples k.
However, in many applications involving rigid-body motion, the measured data corresponds to position and orientation coordinates of the rigid body, instead of translational and rotational velocity vectors.In this case, one option is to use the following two-step approach: (1) estimate the rigidbody velocities from the measured position and orientation coordinates; (2) apply OCP ( 25)-( 26) to the estimated velocity vectors.A downside of this first option is that integration of the optimized velocity vectors to reconstruct position and orientation coordinates will inevitably result in integration drift, such that there will be a deviation between the measured and reconstructed positions and orientations.
A second option is to adapt optimization problem ( 25)-( 26) to include position or orientation measurements.Below, this is worked out assuming the measured positions are given by position vectors and the measured orientations by rotation matrices.The proposed approach can be adapted to deal with other orientation representations such as quaternions.
Including position measurements for case c = v: The reconstructed position vectors p k at each sample are introduced as additional variables in the OCP.The trajectory reconstruction cost ∆c MSE (c, c meas ) in ( 23) is now replaced by the MSE error between the reconstructed position p k and the measured position p meas k over k = 1...N : The reconstructed positions are defined by extending the integrator in (21) with a fourth row and column: The effect of these additions is that c k , which in this case signifies the translational velocity in the moving frame, is used to update the position p k .
Including orientation measurements for case c = ω: The workflow is similar as for position.Variables for the reconstructed rotation matrices R k at each sample are introduced.The trajectory reconstruction cost ∆c MSE (c, c meas ) in ( 23) is replaced by a first-order approximation of the mean-squared rotation angle [41] between the reconstructed rotation matrix R k and the measured rotation matrix R meas k over k = 1...N : where ∥•∥ Fu corresponds to a Frobenius norm that is applied to the upper-triangular part of the matrix.The reconstructed rotation matrices R k are defined by integrating c k , which now signifies the rotational velocity, in an additional reconstruction constraint: Finally, adding the constraint R T 1 R 1 = I 3×3 ensures that orthonormality is preserved throughout the trajectory.

IV. INVARIANT DESCRIPTOR FOR SCREW TRAJECTORIES
This section derives a local invariant descriptor for a screw trajectory s(ξ), parameterized with respect to a progress variable ξ: where a is the direction part of the screw and b the moment part.Again we omit the explicit dependency on ξ for conciseness while the derivative of s with respect to the progress ξ is denoted as: s ′ = ds dξ .

A. Definition of invariants
The screw invariants are based on a generalization of screw theory for rigid-body kinematics and statics.In 1900, Ball [35] published a complete theory of screws based on the theorems of Chasles and Poinsot.Chasles' theorem states that any rigidbody displacement can always be represented as a rotation and translation along an axis in space.Poinsot's theorem is the dual to Chasles' theorem and it states that all forces and moments acting upon a rigid body can always be replaced by a single force and moment along an axis in space.This duality allows us to calculate the screw invariants of a wrench trajectory similar to how it was done for a twist trajectory in [25].Hence, s can represent either a twist t or a wrench w.In screw theory, the screw axis refers to the axis along which the twist or wrench can be expressed [42], but commonly it is also referred to as the instantaneous screw axis (ISA) to highlight its instantaneous nature.
The screw invariants are expressed in a local moving frame, referred to as the Instantaneous Screw Axis frame (ISA frame), which is constructed as follows.The orientation of the ISA frame is completely determined by a, the direction part of the screw, i.e. ω for a twist and f for a wrench.Hence, this orientation, represented by the rotation matrix R, is calculated in a similar way as the FS frame using (15) after replacing c by a. Now, the tangent e t is along the ISA, while the binormal e b is in the direction in which the ISA is instantaneously rotating.On the other hand, the origin of the ISA frame p is defined as the point on the ISA about which the ISA is instantaneously rotating.An analytic procedure to derive this origin was proposed in [25].We summarize it here as a twostep procedure, illustrate it schematically in Fig. 3a, and refer to [25] for more detail.First, the point p ⊥ on the ISA closest to the origin of the reference frame is determined: Second, the parallel displacement p ∥ along the ISA (from p ⊥ to p ) is defined [25]: such that the origin of the ISA frame is found by their addition: As opposed to the FS frame, the ISA frame is a complete moving frame having both a position and orientation, which can be expressed using the pose matrix T : With the ISA frame determined, the screw invariants can be defined.The first two invariants a and b are defined, respectively, as the components of a and b along the ISA: These two invariants are referred to as object invariants since they relate to the object: for motion, they correspond to the rotation and translation of the object along the ISA, while for force, they correspond to the force and moment applied to the object along the ISA.The screw trajectory is reconstructed from the object invariants using a screw transformation matrix S(T ) as in (7): The first-order kinematics T ′ of the ISA frame T is completely determined by four invariants: Here, ω κ and ω τ correspond to the curvature rate and torsion rate of the vector a, which is the direction part of the screw s.
Hence, the invariants ω κ and ω τ are calculated using similar formulas as in (19), with c replaced by a.This is evident because, for all ξ, the orientation part R in (36) is calculated in exactly the same way as R in (15) for the FS frame.Hence, also its derivative is the same.The invariants v b and v t correspond to the translation of the origin of the ISA frame along e b and e t and are defined as [25]: Equation ( 39) can be written compactly using (4): where i is referred to as the moving frame invariant since it models the first-order kinematics of the ISA frame.Fig. 3 visualizes the definition of the moving frame T , the object invariants a and b, and the moving frame invariant i.
Singularities: The singularities of the screw invariants bear a strong relationship with the singularities of the vector invariants.The first type of singularity occurs when the direction vector is zero (a = 0).The object invariant a is zero while the moving frame invariants ω κ , ω τ , v t , and v b are undefined.In other words, the complete moving frame is arbitrary.In this case, the moment vector b, i.e., v for twist and m for wrench, can still be used to construct the orientation of the moving frame using the vector invariants defined in Section III.This singular case is also referred to as infinite pitch or pure moment for motion and force, respectively.
The second type of singularity occurs when the screw axis has an unchanging direction (a × a ′ = 0) so that it remains parallel to itself.As a result, the invariant ω κ is zero, while the other moving frame invariants ω τ , v t , and v b are undefined.The perpendicular distance p ⊥ is still defined.In other words, the direction of the ISA and perpendicular distance to the ISA are well-defined in this case, but the parallel distance p ∥ and normal and binormal axes are arbitrary.In this case, the derivative of the perpendicular distance p ′ ⊥ can be used instead of a to construct the normal and binormal axes.Thus, it is possible to determine the complete orientation and perpendicular distance of the ISA frame, but not its position along the ISA since p ∥ is undefined.
Invariant properties: All the invariant properties of the vector invariants in Section III also hold for the screw invariants.The key advantage of the screw invariants compared to vector invariants is their additional invariance to the choice of the location of the reference point on the object.In other words, screw invariants are invariant to the choices of both the world frame {w} and body frame {b} in which the screw coordinates s may be expressed.This property is also referred to as bi-invariance [36].
Trajectory reconstruction: Reconstructing the trajectory s(ξ) from the invariant descriptor ( s(ξ), i(ξ)) requires numerical integration of the differential equation (41).Discretizing the problem again, the moving frame T k+1 at sample k + 1 can be found from T k at sample k as follows: where we assume that i k remains constant over the integration step ∆ξ.The matrix exponential operator exp(•) maps the matrix [i] × ∆ξ which is part of the Lie algebra se(3) into the corresponding change in pose which is part of the Lie group SE(3), and has a closed-form expression under the given assumptions [36].Given the initial value for T 1 at k = 1, the moving frame T k can be reconstructed for all remaining samples k = 2...N .From the moving frames, the corresponding screw trajectory s k is retrieved using (38): Table III (bottom half) provides an overview of the screw invariants framework with applications to rigid-body twist t and wrench w trajectories.

B. Robust calculation of invariants using optimal control
Calculating screw invariants s and i with the analytical formulas has similar problems related to singularities, sensitivity to noise, and integration drift during reconstruction, as mentioned in Section III-B for vector invariants.
To address these problems, we again formulate an OCP, of which the objective is to find the screw invariants that reconstruct the given measured screw trajectory in a robust way.This is reflected by introducing a cost ∆s k that contains the difference between the reconstructed s k and measured screw trajectory s meas k at sample k: Similarly as in Section III-B, a mean-squared trajectory reconstruction error is formulated, but the errors on the directional part of the screw ∆a MSE and the moment part of the screw ∆b MSE are considered separately: bounded by tolerance ϵ a for the direction part of the screw and tolerance ϵ b for the moment part of the screw.A regularization cost is introduced for the moving frame invariants to achieve a more stable and smooth evolution of the moving frame T in the presence of measurement noise and singularities: where L signifies a chosen length scale with units [m] to properly compare rotational and translational invariants.
The complete OCP to calculate the screw invariants s and i from a given sequence of N measured screws s meas k is: minimize subject to: Besides the trajectory error constraints ( 45)-( 46), the constraints consist of the trajectory reconstruction equations ( 42)-( 43), supplemented with an orthonormality constraint on the initial moving frame T 1 .
Initialization: Similarly as for vector invariants, an invariant initialization of the variables in the OCP can be done using either the analytical formulas or an approach based on an average moving frame [40].The latter is summarized as follows.The orientation of the moving frames T k is initialized using the same approach as for vector trajectories, based on the directional part of the measured screw trajectory s meas k .The origin of the moving frames T k is initialized using the average intersection point of all the instantaneous screw axes corresponding to the measured screw trajectory s meas [40].
The moving frame invariants ω κ , ω τ , v b , and v t are initialized with zero values.The object invariants a k and b k are initialized by transforming the screw trajectory s meas k to the average moving frame and selecting the components along the first axis for each sample k.
Application to motion and force data: The calculation of screw invariants using OCP (48)-(49) can directly be applied to twist trajectories t(ξ) and wrench trajectories w(ξ).
Similarly as in Section III-B, if the measured motion data consists of positions and orientations of the rigid body, instead of twists, the optimization problem (48)-( 49) can be adapted to include the object's pose.Below, this is worked out assuming the positions and orientations are given by homogeneous transformation matrices.The proposed approach can be adapted to other pose representations, such as dual quaternions.
Including pose measurements for case s = t: Additional variables are introduced for the reconstructed pose matrices The trajectory reconstruction costs ∆a MSE and ∆b MSE in ( 45)-( 46) are replaced by the MSE errors between the reconstructed and the measured rotation matrices and positions: The reconstructed pose matrices T k are defined by integrating s k , in this case signifying the object's twist, in additional reconstruction constraints: Finally, adding an orthonormality constraint on the initial pose T 1 , like in (49), ensures that orthonormality is preserved.

V. APPLICATION TO 3D CONTOUR FOLLOWING
This section demonstrates the use and benefits of screw and vector invariant descriptors for both motion and force trajectories by analyzing a human-demonstrated 3D contour following task, shown in Fig. 4. The human operator holds the tool while following the contour and pressing the tool's five contact wheels (diameter 22 mm) against the contour.
The contour consists of the edge between two contact surfaces that form an angle of 90 • in each point along the contour (see Fig. 4c).Each of the five contact wheels takes away one degree of freedom for the motion of the tool relative to the contour, while creating one degree of freedom to apply a contact force.Hence, every point along the contour is characterized by a 1-dof vector space of possible twists and a 5-dof vector space of possible contact wrenches.This makes the application interesting, as we expect to obtain a very good repeatability of the invariant descriptors for motion among different demonstrations, while we expect to see the effect of human variation among the demonstrations in the invariant descriptors for force.Furthermore, the contour is designed to have a symmetric curvature profile and an antisymmetric torsion profile, both with respect to its midpoint.The contour starts and ends with a straight line segment (no curvature or torsion).

A. Experimental set-up
An HTC VIVE motion capture system, consisting of a tracker mounted on the tool and a camera (not shown in Fig. 4), records the pose trajectory of the tool with respect to a reference frame attached to the camera.This frame is chosen as the world frame {w}.The system records the pose trajectory with a frequency of 200 Hz and an accuracy in the order of a few millimeters and a few degrees.The contact force/torque is measured with the same frequency using a 6axis JR3 force/torque sensor.The nominal accuracy of this force/torque sensor is ±1% of its standard measurement range which is ±800 N along the central axis (Z-axis) and ±400 N in the horizontal plane (XY-plane) for the force, and ±24 Nm in all directions for the moment.
Three frames are rigidly attached to the tool: the motion tracker frame {tr}, the force/torque sensor frame {f s}, and the tool-center-point (TCP) frame {tcp}.The {tcp} frame is defined such that, when the tool is tracking a straight-line segment, the origin of {tcp} corresponds to the intersection of the contour's edge and the symmetry axis of the cylindrical part of the tool (see Fig. 4c).One axis of the {tcp} frame is along the edge and the remaining two axes are along the two surfaces.(For completeness: if the edge has curvature and/or torsion, the {tcp} frame will deviate slightly from this contourcentered definition due to the complex contact geometry).
To introduce variation in the measurements, the tracker can be physically attached to the tool at different locations, further discussed in Section V-C.However, the location of the force sensor relative to the tool cannot be changed in the current setup.Instead, and without loss of generality, we can artificially transform the physically measured wrenches to a different, virtual frame on the tool using a screw transformation matrix S as in (9).We can take this one step further and in a similar manner transform the physically measured wrenches to the world frame {w}.This way we can, again without loss of generality, simulate an alternative experimental set-up in which the force sensor is fixed to the world instead of the tool, i.e., representing a situation where the force sensor is placed underneath the support to which the contour is fixed.

B. Objectives
A first objective is to show that the invariant descriptors have a physical interpretation and to obtain evidence that they are calculated correctly using the OCPs.It is however very difficult to calculate the ground truth invariants given the complex contact geometry involving five contact wheels with finite radii and finite mutual distances.Instead, we calculated reference values to compare the invariants with, for the limit case of an infinitely small set of contact wheels and a contact wrench which corresponds to a pure force with constant magnitude ∥f ∥ = 25 N, as shown in Figure 4c.In this limit case, the pose of {tcp} is easily derived from the contour's CAD data.The wrench is then also uniquely defined and constant in any frame attached to the tool.From the reference pose and reference wrench we calculate the corresponding invariant descriptors using the same OCPs as for the humandemonstrated data.
A second objective is to confirm the invariant properties of the vector invariants and the screw invariants, as discussed in Sections III and IV, respectively.Furthermore, we want to point out the practical benefits of these properties in terms of relaxed calibration needs.We also want to show the difference between the invariant descriptors depending on the perspective of the input data: as the motion and wrench trajectories can be defined either with respect to the world or with respect to the moving body, the corresponding moving frame invariants will also model the motion of the moving frame with respect to the world or the moving body, and hence they will be different.
A third objective is to confirm that the new formulation of the OCPs for calculating the invariants, as introduced in Sections III-B and IV-B, allows for an intuitive and effective parameter tuning.
A fourth and final objective is to confirm that the trajectories reconstructed from the calculated invariants are accurate and driftless, even when the trajectories are reconstructed in a different region in space.

C. Experiment design
Twelve demonstrations of the 3D contour following task were recorded.According to the experimental protocol, the operator was asked to repeat the task twelve times by starting from approximately the same location (at the end of the initial straight-line segment) and approximately ending at the same location (at the beginning of the final straight-line segment).The operator was also asked to try to maintain a constant pure contact force as shown in Fig. 4c, because this can be shown to be a very good strategy to ensure that the five wheels remain in contact with the contour.Both instructions contributed to limiting human variation in the task execution.
Each of the twelve demonstrations was recorded with the tracker physically attached to the tool in one of six different locations, shown in Fig. 5 (i.e. two demonstrations per tracker location), resulting in significantly different pose trajectories.Additionally, the physically measured contact wrenches were artificially transformed from the force sensor frame {f s} to the tracker frame {tr} used in the demonstration, simulating a change in location of the force sensor and resulting in significantly different wrench trajectories.
The screw invariants are expected to be invariant for the previous variations, but the vector invariants are not.To verify the repeatability of the vector invariants, it is necessary to pick a well-chosen reference point on the tool to express the position and moment trajectories, the same reference point for all demonstrations.The obvious choice was the origin of {tcp}, both for the position and moment trajectories.Hence, for this specific case, the measured pose and wrench trajectories of the twelve demonstrations were additionally transformed to the {tcp} frame.Note that this requires calibration of the origins of {tr} and {f s} with respect to {tcp}, but this is required solely for the vector invariants, not the screw invariants.
The twelve resulting pose trajectories were further modified by changing the pose of the world reference frame relative to the contour, simulating a change in camera viewpoint.Table IV lists the different artificial transformations that were applied to each of the twelve trials together with their corresponding tracker configuration.Changing the world reference frame is not expected to affect the screw or vector invariants for motion.Obviously, such change also does not affect the wrench trajectories and the resulting wrench invariants since the force sensor is attached to the tool, not to the world.
To simulate the alternative experimental set-up in which the force sensor is fixed to the world, the twelve measured wrench trajectories were artificially transformed to a single frame, fixed to the world underneath the contour, about halfway the contour.The exact location of this frame relative to the contour is not known or calibrated, but this does not affect the screw invariants.
Finally, to show that the trajectories reconstructed from the invariants are accurate and driftless, the motion and force trajectories were reconstructed at a new location and were compared with the measured trajectories, after they were globally transformed to the same location.

D. Data processing
This section explains how the screw and vector invariants are calculated from the measured pose and wrench trajectories and the reference pose and wrench trajectories.
Preprocessing: There are two preprocessing steps: segmentation and re-parameterization.
The purpose of segmentation is to extract the motionin-contact part from a complete demonstration, which also includes an approach-to-contact motion, a retract motion, and two stationary parts while in contact, at the start and the end of the actual contour tracking motion.Segmentation is achieved by using a threshold for the magnitude of the force (5 N) to detect contact, and a threshold for the magnitude of the translational velocity of the TCP (0.05 m/s) to remove the stationary parts of the demonstration.
Re-parameterization involves changing the time-based input trajectories to trajectories that are function of an alternative progress variable.As discussed in Section II-C, for tasks with dominant translation like contour following, the integral of the translational velocity is a natural choice.Accordingly, we chose the progress variable ξ as the integral of the TCP's velocity, because it best reflects the arc length along the contour.
Finally, we rescaled the progress variable to a dimensionless variable ranging from 0 to 1, making it independent of the scale of the contour.
Calculation of the invariants: The force invariants were calculated from the wrench measurements using the original OCPs: ( 25)-( 26) for vector invariants and ( 48)-( 49) for screw invariants.The motion invariants were calculated from the pose measurements using the adapted versions of the OCPs using ( 28)- (31) for the vector invariants and (50)-( 52) for the screw invariants.All OCPs were specified using the nonlinear optimization framework CasADi [43] and solved using a primal-dual Newton method with IPOPT [44].
Initialization of the OCP: For initializing the moving frames, we followed the initialization approaches based on an average moving frame as explained in Sections III-B and IV-B with a minor variation.Due to the symmetry and antisymmetry in the contour, averages calculated over the trajectory might become approximately zero.To avoid confusion in determining the signs, we calculated the average over the last twothirds of the trajectory.The motion invariants are calculated using pose measurements as input.However, the initialization routine uses velocity and twist data as input.Hence, for motion, numerical differentiation is performed to estimate the velocities and twists as inputs for the initialization routine.
Tuning of the OCPs: Vector and screw invariants for motion and force were calculated using measured position, orientation, force, and moment trajectories as inputs.For every type of trajectory, a tolerance ϵ related to the desired accuracy of the reconstructed trajectory had to be chosen, such as for example in (24) for the vector trajectories.To maintain an intuitive relation between the values for these tolerances and the accuracy of the measurement systems, we expressed these tolerances in the respective sensor frames.All reported results were obtained using the following intuitive values: ϵ p = 2 mm for position and ϵ R = 2 • for orientation, which corresponded to a rough estimation of the accuracy of the motion capture system, and ϵ f = 0.8 N for force and ϵ m = 0.16 Nm for moment, which corresponded to the estimated noise levels of the force/torque measurements.To calculate screw invariants, a length scale L must be chosen to properly compare rotational and translational invariants.For the application we chose L = 0.5 m, i.e. the length of the contour's edge.
Reconstructing trajectories from invariants: The pose trajectory T (ξ) of the tool and the contact wrench trajectory w(ξ) were reconstructed from the screw invariant descriptors at a new location in space.This was done by supplying new initial moving frames T (0) for both motion and wrench, and a new initial pose of the tool T (0) in the trajectory reconstruction equations ( 42), ( 43) and (52).The frames were specified by globally transforming the original calculated initial frames to the new location.Time-based trajectories T (t) and w(t) were obtained by re-applying the time profile ξ(t) that was extracted from the demonstration using the inverse of the reparameterization in (11).respect to the world (Fig. 6), force invariants with respect to the tool (Fig. 7), and force invariants with respect to the world (Fig. 8).Additional results for other cases can be generated in the provided software [34].Each plot contains twelve thin lines corresponding to each of the twelve demonstrations.The thick dotted blue line represents the reference values for the limit case.The variation among the demonstrations is given by the two-sigma band depicted in gray.For all invariants, the dimensionless arc length ξ defined in Section V-D is chosen as the progress variable.Figure 10 visualizes the corresponding evolution of the different moving frames for one demonstration.Note that in the subfigures of Fig. 10, the contour is only positioned approximately3 because its true pose was not recorded.Figure 11 visualizes motion and force trajectories reconstructed from the invariants at a new location.Below, the results are discussed in detail.

E. Results and discussion
Motion invariants with respect to world: Fig. 6a shows the screw invariants for motion.The object invariants a and b describe the rotation and translation of the object along the screw axis, while the moving frame invariants ω κ , ω τ , v b , and v t describe the rotation and translation of the moving frame attached to the screw axis.Recalling that each demonstration has different tracker and world frames, the relatively small variations between demonstrations confirm the invariance properties of the screw invariants.Even in the presence of human variations, relatively small variations between the demonstrations and good correspondence with the reference values from the limit case are obtained.This is explained by the constrained motion, allowing variations in only one dof.Due to the special design of the contour, clear symmetric or antisymmetric profiles with respect to the midpoint were expected and were obtained.The motion of the moving frames for one trial is visualized in Figure 10a.
Fig. 6b shows the vector invariants for orientation.According to the analytical formulas, they should correspond to the invariants for the directional component of the screw for motion (top row in Fig. 6a).The reason why they are not exactly the same is explained by the OCP-based calculation approach.In this approach, the calculation for orientation and translation are coupled into a single OCP for the screw invariants and hence they are subject to a joint regularization.Nevertheless, the overall evolution of the vector invariants for orientation and the invariants for the directional component of the screw are still similar.
Turning to the vector invariants for translation, Fig. 6c confirms that they are not invariant for changes of the tracker location, serving as the reference point to define the translation of the tool.This non-invariance can be resolved by picking the same reference point for all demonstrations, in this case the origin of {tcp}, as shown in Fig. 6d.In this figure, c represents the magnitude of the translational velocity vector, while ω κ and ω τ represent the dimensionless curvature and torsion of the point curve tracked by the origin of {tcp}.From all the figures, Fig. 6d exhibits the best correspondence between the reference limit case and the demonstrations.This is to be expected since the limit case models the translation of the {tcp}-frame along the edge.For this limit case, c is approximately constant (about 0.5 m/-).Its integral (also close to 0.5 m) represents the arc length of the point curve tracked by the origin of {tcp}.The small values for the torsion rate indicate that this point curve lies, at least locally, approximately in a plane.The small differences between the limit case and the demonstrations are due to the finite dimensions of the tracking tool and due to measurement errors.The moving FS frames corresponding to these vector invariants for one demonstration are visualized in Fig. 10d.
Given the geometry-based definition of the invariants, they can be used to segment the trajectory into meaningful parts.For example, in the screw invariants in Fig. 6a, we notice two inflection points at ξ ≈ 0.2 and at ξ ≈ 0.8: while the rotational velocity a changes sign at these progress values, the peaks in v b show that the position of the ISA is rapidly shifting from a center of curvature on one side of the contour to the other side of the contour.This is confirmed by the visualization of the corresponding ISAs in Fig. 10a.
Force invariants with respect to tool: More variation (a) Screw invariants calculated from wrench transformed to tracker frame (b) Vector invariants calculated from force transformed to tracker frame ) Vector invariants calculated from moment transformed to tracker frame ) Vector invariants calculated from moment transformed to tool center frame Fig. 7. Screw and vector invariants calculated from measured force and moment trajectories, transformed to either the tracker or the tool center frame.
among demonstrations was expected in the wrench invariants because, although the operator was instructed to apply a pure force in a specified direction and through a specific point of the tool, the operator could at any time apply a wrench in a local 5-dof space.The higher noise levels associated with force measurements may cause additional variation compared to the motion invariants.Note that a pure force that remains parallel to itself corresponds to the second type of singularity, mentioned in Section IV.This singularity applies to both vector and screw invariants.
Fig. 7a shows the screw invariants.Recalling that each demonstration has a different transformed force sensor frame, these plots confirm that the screw invariants are invariant for such transformation.The plot of invariant a shows that the operator was able to maintain relatively constant magnitudes of the force along the contour, although the magnitudes varied between demonstrations.The plots of ω κ and ω τ reveal a limited change of direction of the force (hence of the ISA), almost in a single plane (ω τ is very small).This plane is orthogonal to the edge of the contour.This can be found  by analyzing the direction of the ISAs and is confirmed in Fig. 10b, where the first (red) and second axes (green) of the successive moving frames lie approximately in a plane that is orthogonal to the tangent of the contour.
The second row of Fig. 7a shows that there is no significant moment about the ISA (b is very small) and that the origin of the moving frame does not translate much, both parallel to and along the ISA (both v b and v t remain very small).The latter can also be seen in Fig. 10b.Fig. 7b shows the vector invariants for force, which are nearly identical to the directional component of the screw invariants (top row in Fig. 7a).As for the vector invariants for moment, Fig. 7c confirms that they are not invariant for a change in position of the force sensor frame, resulting in moment magnitudes up to 2.5 Nm.Similarly as for translational velocity, this can be resolved by picking the TCP as the common reference point for expressing the moments in all demonstrations, as shown in Fig. 7d.Since these moment values are relatively close to zero w.r.t. the specified tolerance ϵ m = 0.16 Nm, they can be perfectly described by a stationary Frenet-Serret frame.Therefore, the solutions of ω κ and ω τ are below 10 −8 , the stopping criterion in the OCP solver.
Force invariants with respect to world: Figure 8 shows the screw invariants for the simulated, alternative experimental set-up in which the force sensor is assumed fixed to the world.This simulation, in particular the transformation of the physical force measurements to the world, is subject to two additional errors: the calibration between the force sensor and the motion tracker, and the camera measurements.Again, the invariants are shown to be insensitive to this transformation.Fig. 10c shows the evolution of the corresponding moving frame.During execution, the force ISA remains approximately perpendicular to the contour edge while it approximately intersects the contour edge.To do this for the entire trajectory, the force axis has to translate along the contour with a velocity v b , which is approximately equal to the translation of the TCP per dimensionless arc length (about 0.5 m/-).Integrating this value over the total horizon again results in the length of the entire contour (∼0.5 m).
To further study the variation in the force invariants with respect to the world, we replaced the real wrench measurements in each demonstration with a constant simulated wrench w = 25 0 cos( π 4 ) sin( π 4 ) 0 0 0 T expressed in the TCP frame.By using this simulated wrench, the effects of human variations and wrench measurement noise were eliminated.Consequently, the corresponding force invariants with respect to the world are only influenced by errors in the calibration between the TCP and the motion tracker, and noise in the motion measurements.Invariant descriptors of this case are shown in Fig. 9.We observe greater repeatability in force invariants compared to Fig. 8 as a result of eliminating the human variation and wrench measurement noise.This analysis confirms that the larger variation in the force invariants compared to motion invariants is mainly due to human variation in a 5-DOF space, not due to errors in the computation.
Reconstructed motion and force trajectories:   trials.Table V compares these reconstructed trajectories with the measured trajectories on which a global transformation has been applied so that they are in the same location.The RMS-differences between the trajectories were found to be in agreement with the chosen tolerances for the trajectory accuracy constraint in the OCPs.The small difference is due to the convergence tolerance of the numerical solver of the OCPs, which was set to 10 −8 .These results confirm that the trajectories reconstructed from the calculated invariants have the accuracy specified in the OCPs and are also driftless.Evidently, these results also hold for the special case of a reconstruction at the same location as the measurements.

VI. APPLICATION TO PEG-ON-HOLE ALIGNMENT
This section examines motion and force trajectories resulting from a human-demonstrated peg-on-hole task in which the operator holds a tool to which the peg with diameter 50 mm is attached, as shown in Fig. 12.The task aims to ensure that peg and hole are aligned before peg insertion starts, even if the hole is inaccurately located.To accomplish this, the operator first establishes a stable three-point contact between the hole and a purposefully misaligned peg, as shown in Fig. 12c.This contact involves two rim-rim contact points and one surfacerim contact point.We recorded 12 demonstrations of the alignment, each starting from such initial three-point contact and finishing at full alignment.The operator was instructed to perform alignment while maintaining the three-point contact.All demonstrations were performed by the same operator, but in two batches, respectively of 7 and of 5 trials, at different points in time.Although this task appears to involve a pure rotation, previous theoretical studies have shown that, apart from a rotation, the motion involves a small, but non-zero translation, see e.g.[46].It is evident that, due to the three-point contact and neglecting friction, this task involves two instantaneous 3dof vector spaces: one for the wrench and one for the motion.So, the operator can use different motion and wrench strategies during the demonstration.

A. Experimental set-up and data processing
The motion capture system is the same as for the contour following task of Section V.The force/torque sensor is again a 6-axis JR3 sensor with a nominal accuracy of ±1%, but its measurement ranges are smaller: ±200 N along Z and ±100 N along X and Y for force, and ±5 Nm in all directions for moment.The same three frames as on the contour following tool were defined on the demonstration tool here (see Fig. 12b), with the TCP frame located at the tip of the peg.
Unlike in the contour following application, we did not use different tracker positions on the tool and did not artificially transform the pose measurements by changing the world frame, because invariance with respect to such variations was sufficiently proven in the contour following experiment.We also did not transform the wrench measurements to different locations on the tool, but we did transform them to a frame fixed to the world to investigate the wrench invariants with respect to the world.
Data was processed in roughly the same way as in the contour following application, with the following changes or additions.The segmentation thresholds were set to 1 N and 0.35 rad/s for force and angular velocity, respectively.Since the task is predominantly rotational, the integral of the rotational velocity was chosen as the progress variable to reparameterize the time-based input trajectories, and it was also rescaled to a dimensionless variable ranging from 0 to 1.The initialization of the OCPs based on average moving frames (Section III-B and IV-B) was done on the complete measured trajectories.Furthermore, since the force/torque sensor was more accurate than the one used for contour following (same relative accuracy for a smaller range), we lowered tolerances ϵ f and ϵ m in the OCPs to 0.3 N and 0.1 Nm, respectively.
Finally, to investigate if the operator used different wrench strategies in the two batches, we calculated the invariants corresponding to the average of the wrench trajectories for each of the two batches and compared them to the invariants corresponding to the average of all demonstrations.This was ) Screw invariants for force with respect to world (c) First two screw invariants for the averaged wrench trajectories of measurement batch 1 and 2, with respect to either world or tool.done both with respect to the tool and with respect to the world.Three remarks have to made here: 1) we took the average of the wrench trajectories rather than the average of the invariants, because averaging is a linear operation, but the subsequent calculation of the invariants is not; 2) we can just average the measured wrench trajectories, because all wrench trajectories were recorded at the same physical sensor location on the tool and were additionally transformed to the same virtual location in the world 4 ; 3) because averaged wrench trajectories are less noisy, we could further lower ϵ f to 0.17 N and ϵ m to 0.05 Nm.

B. Results and discussion
Figure 13 and 14 summarize the experimental results, focusing on screw invariants for motion and force.Additional results for the vector invariants can be generated in the provided software [34].
Motion invariants with respect to world: Fig. 13a shows the screw invariants of motion.The object invariant a corresponds to the relative rotation of the peg with respect to the hole around the screw axis.Its value is almost constant because the re-parameterization is based on the rotation and it lies around 1 rad/-.This makes sense since the integral of a corresponds to the initial misalignment angle, which is about 1 rad.Object invariant b is small, meaning there is almost no translation along the screw axis.The moving frame invariants ω κ , ω τ , v b , and v t are zero or near zero, indicating that the moving frame is stationary.Hence, with the chosen values of the tolerances in the OCP, the relative motion between peg and hole can be explained by a pure rotation about a fixed axis.As a result, the small theoretical translation mentioned earlier vanishes due to the chosen tolerances.This conclusion holds for all twelve trials, indicating that the operator used the same motion strategy in all trials: a pure alignment of peg and hole achieved by a planar motion, without exciting the two other relative degrees of freedom, being relative rotations about peg or hole axes.This is confirmed by Fig. 14a where all the instantaneous screw axes and frames coincide.
Force invariants with respect to world: Fig. 13b shows the screw invariants for the individual demonstration trials.Object invariant a, which corresponds to the magnitude of the force vector, is remarkably small, going from less than 4N to zero during the alignment.As for the moving frame invariants, only ω κ is significantly different from zero, at least for six of the twelve trials.This indicates that, for those six trials the instantaneous screw axis (i.e. the force vector) is rotating with respect to the world in a plane about a fixed point; for the other six trials the instantaneous screw axis remains fixed to the world, i.e. the force vector does not rotate.This in turn indicates that the operator used different wrench strategies in the trials.
Force invariants with respect to world and tool for averaged wrench trajectories: Based on the previous result, we hypothesized that the operator had used different wrench strategies in the two batches of demonstration trials, because they were executed at different points in time.Figure 13c therefore shows the screw invariants for the averaged wrench trajectories of batch 1 (red), batch 2 (blue), and the whole set of trials (green).We show both the invariants with respect to the world (solid lines) and to the tool (dashed lines).Only the first two rotational invariants are shown, as invariants are zero (v b , v t ) or almost zero (ω τ , b), similarly as in Fig. 13b.Visually, there is an evident distinction between the trials of batch 1 and batch 2: while for batch 1 the average of the trials shows a force vector that rotates with respect to the world and remains fixed with respect to the peg, the opposite is true for the average of the trials of batch 2. To illustrate this further, we plot the resulting moving frames for one trial from batch 1 and one trial from batch 2 in Figures 14b and  14c.These figures clearly show that the ISA is rotating for the trial from batch 1 while the ISA is stationary for the trial from batch 2. A more in-depth study and statistical analysis of these wrench strategies go beyond the scope of this paper, but we conclude that different strategies can lead to successful alignments.This is a valuable result for designing robot force control strategies.It is also interesting to see that it could be obtained in an experiment with only small contact forces.

VII. DISCUSSION AND CONCLUSION
The aim of this work was to introduce local invariant descriptors for force and moment trajectories.Inspired by the invariant descriptors for motion, we exploited the duality between motion and force.We extended existing concepts for motion trajectories to general vector and screw trajectories, and then applied them to force and moment trajectories.
Two types of descriptors were introduced: vector and screw invariants.Vector invariants describe 3D vector trajectories (translational velocities, rotational velocities, forces, and moments) while screw invariants describe 6D screw trajectories (twists and wrenches).Both types of descriptors are invariant for the choice of world reference frame, but the screw invariants are also independent of the reference point on the moving body that is chosen to represent the translational velocity or the moment.Hence, for obtaining the screw invariants there is no need to choose a reference point on the moving body.This is particularly useful in cases where there is no natural choice for such reference point, the reference point can change between executions, or the calibration of the reference point is time consuming.The newly proposed descriptors can be applied to both motion and force trajectories, as opposed to existing invariant descriptors, which were mainly intended for point curves and/or motion trajectories [16], [19]- [27].
Invariant properties allow to recognize and model tasks that were demonstrated in different environments, without the need for accurate calibration.For trajectory reconstruction and generation however, calibration efforts cannot be avoided since for these applications we need initial values as pointed out in the related paragraphs in Sections III and IV.More in detail, we need initial values for the orientation R 1 or location T 1 of the moving frame, and, if we reconstruct the object trajectory up to position and orientation, we also need respective initial values p 1 and R 1 .In practical applications, these initial values are to be derived from environment models or from sensor measurements, either directly or using modelbased state estimators.
Following a similar line of reasoning, the invariants themselves cannot be used to describe or check relations between motion and force, such as reciprocity between a screw twist and a screw wrench, unless the transformation between the respective moving frames in which the motion and force invariants are represented is known.In general, this requires calibration.For example, to check the reciprocity condition, i.e., the condition that no power is dissipated in the contact: f • v + m • ω = 0, it is required that both twist and wrench, t and w, are expressed in the same reference frame and have the same reference point.Therefore, we cannot just plug in the a and b object invariants of the respective screw invariants, because they are defined relative to their respective moving frames, and, in turn, to the frames in which the twist and wrench trajectories were measured.Hence, to study relations between motion and force, we need a calibration between the respective measurement frames.
The vector and screw invariants can be made independent of time, and hence independent of the applied time-based motion profile, by expressing them as a function of a suitable geometric progress variable.Again, this additional invariance is useful for the recognition of a task demonstrated in different situations and for the generalization of task models to new situations, as shown in [15] for the case of motion in free space.However, the way we chose the progress variable in the 3D contour following application introduced an important limitation.We did not follow the approach to choose an invariant progress variable that is independent of the reference point for translation, as suggested in [25] and briefly discussed in Section II-C.Instead, we chose the progress variable as the integral of the translational velocity of the origin of {tcp}.So, based on our understanding of the contact geometry, we handpicked this point as a reference point to define arc length as a representation of task progress, and used it to represent all invariants.The benefit of this approach is that we obtained vector and screw invariants for both motion and force that are interpretable and can be compared with each other and with the reference values of the limit case (Figures 6-8).Definitely, this progress value is a very good choice and probably even the best one for this application, but it would be better to derive a suitable progress value from the measured trajectories in an invariant data-driven way.This is subject of our future work.
Besides analytical procedures to derive the invariant descriptors, an approach was presented based on optimal control, where the main idea is to find the invariant descriptor that reconstructs the given window of measurements.This provides a way to deal with singularities and measurement noise and hence enables a much more robust calculation of the descriptors.Compared to our previous work on the calculation of motion descriptors [33], we improved the stability and repeatability of the calculation by reformulating the OCP.In the reformulated problem, trajectory reconstruction errors became inequality constraints that were bounded by a desired tolerance value.This resulted in more intuitive tuning, as we showed that these tolerance values could be directly related to the accuracy or noise characteristics of the used motion and force sensor.The reformulated OCP also resulted in a more stable behavior of the moving frame.For example, when the values of the noisy moments are below the specified tolerance, they can be described by a stationary moving frame (see Fig. 7d and Fig. 10b).Finally, the reformulation resulted in a smaller number of tuning parameters, further reducing the tuning effort.Overall, we obtained a good repeatability of the invariants, as can be seen in Figures 6-8, even for noisy data.
Continuity of the moving frames is ensured by minimizing the values of the moving frame invariants in the OCPs.This also provides some degree of continuity for the reconstructed trajectories.To enforce a more specific continuity on the reconstructed trajectory, it is advised to include the time profile as variables in the OCPs in order to enforce acceleration and/or jerk limit constraints.These constraints can be either on the spatial trajectory itself or on the robot joint trajectories in case a kinematic model of the robot is available, similar to [15].
The paper contains a practical experiment involving a human-demonstrated 3D contour following task which showcases the use, verifies the invariant properties, and validates the robust calculation of screw and vector invariant descriptors for both motion and force.As expected, this experiment yielded highly repeatable motion descriptors due to the highly constrained motion (in a 1-dof instantaneous twist space), and less repeatable force descriptors due to human variation (in a 5dof instantaneous wrench space).Evidently, apart from human variation, other disturbances were present in this experiment such as sensor inaccuracy and measurement noise.A second experiment, involving human demonstration of a peg-on-hole alignment task, showed how motion and force invariants can be used to detect and characterize different human strategies in contact tasks.
The concepts and algorithms presented in this paper are intended as an initial toolbox [34] to support future work on invariant representations of motion and force trajectories and their use in physics-informed machine learning of forcecontrolled manipulation skills.

APPENDIX RELATION BETWEEN THE KINEMATICS OF THE FS FRAME
AND THE FRENET-SERRET EQUATION In differential geometry, the Frenet-Serret equation is usually written as follows: with ω κ the curvature, ω τ the torsion, and T , N , B the unit vectors of the FS frame.This equation can be shown to be equal to the kinematics of the FS frame as defined in (18): Since the transpose of a skew-symmetric matrix is equal to the negative of the matrix, (55) corresponds to the Frenet-Serret equation in (53).

which is part of the SE( 3 )Fig. 1 .
Fig. 1.Visualization of screw twist t = (ω T v T ) T and screw wrench w = (f T m T ) T in the world frame {w} and in the body frame {b}.

Fig. 2 .
Fig.2.Invariant descriptor for a vector trajectory: (a) object invariant c defined in the moving FS frame given by the unit vectors et, en and e b ; (b) moving frame invariants ωτ and ωκ.The figure depicts the special case where c corresponds to a translational velocity vector with respect to reference frame {w} and where, for visualization purposes, the moving frame is depicted with its origin chosen at the corresponding position along the point trajectory (indicated in blue).

Fig. 3 .
Fig. 3. Invariant descriptor for a screw trajectory: (a) object invariants a and b defined in the moving ISA frame given by the unit vectors et, en, e b , and the origin p = p ⊥ + p ∥ et; (b) the moving frame invariants ωκ, ωτ , v b , and vt.The figure shows the special case where the screw s corresponds to the screw twist t, representing the motion of the rigid body {b} w.r.t {w}.
c = ω for rot.velocity c = v for transl.velocity c = f for force c = m for moment derivatives c ′ and c ′′ orientation frame R = [et en e b ], object trajectory equation c = R c object orientation R equation: R ′ = [c] × R, when c = ω object position p equation: p ′ = c, when c = v

Fig. 4 .
Fig. 4. The 3D contour following setup: (a) demonstration setup showing tool and contour, (b) 3D drawing showing the three assigned frames on the tool: tracker {tr}, force/torque sensor {f s}, and tool center point (TCP) {tcp}, (c) cross-section of the contour showing the reference contact force f .

Fig. 5 .
Fig. 5.The six different configurations of the motion tracker on the tool.

Fig. 8 .
Fig. 8. Screw invariants calculated from measured wrench transformed to the world frame.

Fig. 9 .
Fig. 9. Screw invariants calculated from a simulated wrench in the world frame to study the reduction in variation with respect to Figure 8.

Fig. 10 .
Fig. 10.Moving frames calculated from data of Trial 5: (a) ISA frames calculated from the tool's motion in the world.The red arrows correspond to the ISA itself, which is the first axis of the moving frame.The green and blue arrows correspond to the second and third axes, respectively.The two blue circles indicate the translation of the ISA frames towards the other side of the contour at the inflection points; (b) ISA frames calculated from interaction wrench, relative to the tool; (c) ISA frames calculated from wrench w.r.t the world; (d) Frenet-Serret (FS) frames calculated from the position of the tool center point (TCP).The red arrows correspond to the tangent along the contour, which is the first axis of the moving frame.
Figure 11   visualizes a reconstruction of the motion and force trajectory at a new location starting from the invariants of one of the

Fig. 11 .
Fig. 11.Reconstruction of the motion and wrench trajectory of Trial 7 at a new location given the screw invariants of the measured trajectory seen from the world.The motion of the follower is shown using the position of the TCP and the orientation of the follower.The force on the follower is depicted by the red arrows at the corresponding position of the TCP.The initial ISA frames for motion and wrench, necessary for reconstruction, are also depicted.

Fig. 12 .
Fig. 12. Peg-on-hole setup with defined frames and illustration of the contact.

Fig. 13 .
Fig.13.Screw invariants for motion and force for the peg-on-hole task.

12 Fig. 14 .
Fig. 14.Moving frames for trials in the peg-on-hole alignment task, all seen from the world frame.The ISAs are depicted using red arrows: (a) Moving frames calculated from the tool's motion for Trial 1 (batch 1), showing a stationary ISA; (b) Moving frames calculated from the wrench for Trial 1 (batch 1), showing a moving ISA; (c) Moving frames of the wrench for Trial 12 (batch 2), showing a stationary ISA.

0 −ω κ 0 ω κ 0 −ω τ 0 ω τ 0 
expression, first expand R into its individual columns using the unit vectors defined in (15): e ′ t e ′ n e ′ b = e t e n e b   .(54)Applying the transpose to both sides of (54) results in:

TABLE II OVERVIEW
OF NOTATION USED THROUGHOUT THE PAPER.(IN THIS TABLE a AND b ARE PLACEHOLDERS FOR TWO FRAMES OF INTEREST.)

TABLE III OVERVIEW
OF EQUATIONS FOR VECTOR INVARIANTS (TOP HALF) AND SCREW INVARIANTS (BOTTOM HALF).THE FIRST COLUMN SHOWS THE TYPE OF VECTOR AND SCREW TRAJECTORIES FROM WHICH INVARIANTS CAN BE CALCULATED.IN THE SECOND COLUMN, THE LOCAL MOVING FRAME IS DEFINED IN WHICH THE INVARIANTS ARE EXPRESSED.THE THIRD COLUMN DEFINES THE MOVING FRAME INVARIANT THAT CHARACTERIZES THE MOTION OF THE MOVING FRAME, AND THE OBJECT INVARIANT THAT CHARACTERIZES THE OBJECT'S TRAJECTORY IN THE MOVING FRAME.THE FOURTH COLUMN DETAILS HOW TO RECONSTRUCT THE TRAJECTORY FROM THE INVARIANT DESCRIPTOR, POSSIBLY INCLUDING THE OBJECT'SPOSITION AND ROTATION.

TABLE IV EACH
OF THE TWELVE DEMONSTRATION TRIALS IS LISTED WITH ONE OF THE SIX CONFIGURATIONS OF THE MOTION TRACKER (SEE FIG. 5) AND WITH THE ARTIFICIAL TRANSFORMATION THAT WAS APPLIED TO THE MEASURED POSES TO SIMULATE A LARGE VARIATION IN THE WORLD (I.E.CAMERA) FRAME.THESE TRANSFORMATIONS ARE SPECIFIED USING XYZ EULER ANGLES AND A TRANSLATION VECTOR.

TABLE V RECONSTRUCTION
ERROR FOR THE MOTION AND FORCE TRAJECTORY OF FIGURE 11.THE CHOSEN TOLERANCES ARE GIVEN AS A REFERENCE.