Haptic rendering and interactive simulation using passive midpoint integration

Abstract

We propose a novel framework of haptic rendering and interactive simulation, which, by exploiting midpoint time integration, that is known for its superior energy-conserving property but has not yet been adopted for haptics and interactive simulation, can enforce discrete-time passivity of the simulation effectively in practice, while retaining real-time interactivity due to its being non-iterative. We derive this passive midpoint integration (PMI) simulation for mechanical systems both in maximal coordinates (i.e. in SE(3)) and in generalized coordinates (i.e. in ${\Re }^n$ ), with some potential actions as well to implement joint articulation, constraints, compliance, and so on. We also fully incorporate multi-point Coulomb frictional contact into them via the PMI-LCP (linear complementarity problem) formulation. The proposed PMI-based simulation framework is applied to some illustrative examples to demonstrate its advantages: (1) haptic rendering of a peg-in-hole task, where very light/stiff articulated objects can be simulated with multi-point contact; (2) haptic interaction with a flexible beam, where marginally stable/lossless behavior (i.e. vibration) can be stably emulated; and (3) under-actuated tendon-driven hand grasping, where mixed maximal-generalized coordinates are used with very light/stiff fingers.

1. Introduction

With many research groups and companies around the world aiming their efforts at blurring the boundary between the real world and virtual worlds in the form of, for example, virtual reality, augmented reality, mixed reality, and so on, providing realistic and believable artificial visual, auditory, olfactory, and mechanical experiences to users in real time is now increasingly in demand. Among these modalities, the artificial mechanical experience, or haptic experience, possesses its own peculiar challenges. Specifically, as compared to other virtual sensory modalities, it requires a much faster update rate while being under more stringent requirements for stability within the virtual world or with the human users.

Numerous strong results and methodologies have been proposed for realizing this virtual mechanical experience by different research communities. On the one hand, in the field of haptics, it is pursued under the name haptic rendering, for which some avatar of the haptic device (e.g. god-object (Zilles and Salisbury, 1995), virtual proxy (Colgate et al., 1995; Lee et al., 2012)) and the virtual objects are simulated together, typically at a rather fast speed (e.g. 500 Hz to 1 kHz), with their complexity encompassing the simplest one-dimensional virtual wall (e.g. Colgate and Schenkel, 1997; Lee and Huang, 2008) to complex deformable objects such as human organs (e.g. Delingette, 1998; Duriez et al., 2006; Garre and Otaduy, 2010). On the other hand, in the field of computer animation (Baraff, 1996; Garre et al., 2011;Müller et al., 2005; Weinstein et al., 2006), it is often called by the name interactive simulation, which typically deals with virtual objects and entities more complex than those of haptic rendering, with a slower update rate (e.g. 10 Hz to 50 Hz), which is still fast enough for users to interact in real time. By utilizing a multi-rate interface (e.g. Duriez et al., 2006; Fotoohi et al., 2007; Lee and Huang, 2010), it is also possible to combine the interactive simulation with haptic rendering.

The core of any frameworks in haptic rendering and interactive simulation of dynamical systems ^* is the time integration, and the majority of state-of-the-art frameworks adopt symplectic Euler integration (SEI (Hairer et al., 2006), also often referred to as the semi-implicit Euler method in the field of haptics and computer animation: see e.g., Cline and Pai, 2003; Müller et al., 2005; Weinstein et al., 2006) or implicit Euler integration (IEI: e.g., Courtecuisse et al., 2015; Duriez et al., 2006; Otaduy and Lin, 2006; Otaduy et al., 2009). Widely used commercial physics engines and interactive simulation tools also rely on SEI (e.g. CHAI-3D ¹ based on ODE, ² and VREP ³ based on ODE or Bullet ⁴ ). Implicit Euler integration, in its exact form, is known to be (strictly) stable (e.g. Chung and Hulbert, 1993; Hairer et al., 2006), allowing for a wide range of parameters (e.g. lighter mass with stiffer spring for transparency) and a longer integration interval (i.e. faster simulation), yet is typically implemented only approximately for interactively-fast running (e.g. linearization (Duriez et al., 2006; Otaduy and Lin, 2006), explicitization (Krysl, 2005, 2006)). This IEI is also in general strictly dissipative, and thus cannot simulate lossless or energy-conserving behavior (e.g. vibrating flexible objects). The SEI is usually simpler to implement (and thus faster) than the IEI due to its semi-implicit nature. This SEI, however, is known to be only conditionally stable (e.g. Hairer et al., 2006), which is fairly detrimental in practice, as one cannot freely tune and explore parameters of virtual objects to achieve (possibly changing) target performance independently from other parts of the virtual world.

On the other hand, passivity has been one of the most powerful tools in the fields of haptics, telerobotics and interactive robotics for enforcing stability among dynamic_ entities coupled via feedback interconnection between, for example, robots and humans, robots and robots, robots and the environment (Buerger and Hogan, 2007; Colgate and Schenkel, 1997; Lee and Huang, 2010; Lee and Li, 2005; Lee and Spong, 2006). This passivity is particularly desirable for haptic rendering and interactive simulation, since, as pointed out in Brown and Colgate (1998), it allows for (1) stable simulation with any parameters and update rate (e.g. arbitrarily light/stiff simulation with coarse integration intervals), and (2) modular construction of provably stable virtual worlds by simply combining individual virtual objects, each designed to be passive. Information delays among virtual objects or virtual worlds can also be addressed in a straightforward manner using this passivity formalism (e.g. for multi-user haptics; Huang and Lee, 2013; Rakhsha and Constantinescu, 2014).

In this paper, we propose a novel passivity-based framework of haptic rendering and interactive simulation. The key differentiating and innovative aspect of our framework is our adoption of the midpoint time integration, whose energy-conservation ability is well known in the field of structural mechanics (e.g. Bathe, 2007; Chung and Hulbert, 1993; Krysl, 2005, 2006; Kuhl and Crisfield, 1999; Simo and Wong, 1991; Simo et al., 1992), but which, to the best of our knowledge, has not been adopted so far for haptic rendering and interactive simulation. This missing of the midpoint integration is particularly surprising for the case of haptic rendering, as its energy-conserving property naturally suggests the possibility of discrete-time (lossless) passive simulation, which was proposed as one of the ultimate goals of haptic rendering in Brown and Colgate (1997, 1998), yet, has been dodging so far any attempts to attain that. Of course, this energy-conserving property is only for a single system, and the passivity property, which is concerned with their interconnection with other objects via power-conjugate input–output mapping, has not been studied in the field of structural mechanics. It was however independently shown in Lee and Huang (2008) and Lee et al. (2012) that, by directly enforcing passivity in the discrete-time domain in a manner similar to that for the derivation of variational integrators (e.g. Marsden and West, 2001), it is indeed possible to develop a passive integrator, which is called the non-iterative passive mechanical integrator (NPMI) and which turns out to be a midpoint integrator. This NPMI was derived in Lee and Huang (2008) and Lee et al. (2012), however, it is limited to only mechanical systems which are linear or written in generalized coordinates with no intermittent and unilateral contact with other objects. Recently, we extended this NPMI for mechanical systems in maximal coordinates with a linear complementarity problem (LCP) formulation of Coulomb friction contact (Kim et al., 2017).

In this paper, we unify and generalize these NPMI-related results by Lee and Huang (2008), Lee et al. (2012), and Kim et al. (2017) as a midpoint integration framework which can deal with linear and nonlinear mechanical systems expressed both in generalized and maximal coordinates under unilateral/intermittent multi-point Coulomb friction contacts formulated via LCP, all the while effectively enforcing discrete-time passivity of the simulation and also its real-time interactivity by preserving the non-iterative nature of the NPMI. Some approximations, inevitable in rendering the algorithms interactively fast, are also adopted here, as in most other haptic rendering and interactive simulation results (e.g. Cline and Pai, 2003; Courtecuisse et al., 2015; Duriez et al., 2006; Krysl, 2005, 2006; Müller et al., 2005; Otaduy and Lin, 2006; Otaduy et al., 2009; Weinstein et al., 2006). Even so, our proposed algorithms can still exhibit fairly passive behavior, as demonstrated by our validating simulation examples (see Section 5) and also as similarly concluded in Krysl (2005, 2006) for midpoint integration with some approximations. For this reason, we call our integrator a passive midpoint integrator (PMI) and our framework, the PMI-based haptic rendering and interactive simulation.

Due to its passivity enforcement and implementation flexibility, our proposed PMI-based framework possesses the following advantages.

1. It can stably simulate mechanical objects with a wide range of parameters (e.g. any positive-definite mass, stiffness and damping gains, and varying update rate possible), thereby allowing for easy and modular tuning of virtual objects while also providing the capacity to adjust the discrete-time damping dissipation necessary to render the simulation more or less stable depending on the target performance and behavior.

2. Consequently, it can simulate very light/stiff objects, which is desired for sharper/transparent haptic rendering (see the peg-in-hole simulation in Section 5.1) or necessary to simulate some classes of mechanical systems (e.g. an under-actuated tendon-driven (UATD) hand with light fingers and stiff tendons; see Section 5.3), but typically cannot be stably rendered by (conditionally stable) SEI-based frameworks (e.g. CHAI-3D).

3. It can simulate marginally stable/energy-conserving behavior such as the vibration of flexible objects (see Section 5.2 and also Section 2.3) or elastic contact between two rigid objects with no extra restitution models (see Section 3.3), which typically cannot be achieved by dissipative/decaying IEI-based frameworks, with its implementation complexity at least comparable to that of IEI.

4. It allows us to use both the maximal and generalized coordinates and even their combination, substantially streamlining the simulation of some mechanical systems such as the UATD hand composed of fingers with tendon routing and a rigid body palm (see Section 5.3) or the flexible beam with a rigid plate connected to haptic devices (see Section 5.2), which, so far, are not so readily integrated in CHAI-3D or VREP.

To our knowledge, this paper provides a solution for the first time to the long-standing problem posed in Brown and Colgate (1997, 1998) of how to attain passivity and interactivity of simulation at the same time.

The rest of the paper is organized as follows. Section 2 contains some preliminary materials, including the definition of passivity, the kinematics and dynamics of rigid bodies in SE(3) and robots, and a summary of the NPMI development for linear mechanical systems in Lee and Huang (2008). We then derive the PMI-based simulation framework in maximal coordinates in Section 3 and in the generalized coordinates in Section 4 with the PMI-LCP formulation of Coulomb friction point contact modeling for each setting. Illustrative examples, namely a peg-in-hole task, flexible beam haptic interaction and UATD hand grasping, are presented in Section 5 to validate and demonstrate our PMI-based framework. A summary and some comments on future research directions are then given in Section 6.

2. Preliminary materials

2.1. Mechanical systems and passivity

The virtual worlds we aim to construct in this paper are mainly composed of multiple articulated rigid bodies (Weinstein et al., 2006) with some constraints (e.g. no penetration, rotation/pivot joints, etc.), elastic elements (e.g. springs), viscous damping elements, and Coulomb friction between them. A part of these mechanical objects is typically interfaced with some human commanding device (e.g. a haptic device), which will not be discussed here in detail; instead, we refer readers to prior results on the (possibly multi-rate) virtual coupling applicable to our PMI simulation (e.g. Huang and Lee, 2011; Lee and Huang, 2010; Lee et al., 2012) and also to Section 5 for some application examples. Large-scale deformable objects (e.g. human organs; Delingette, 1998) are not the main concern of this paper either and we leave them and their salient issues (e.g. model reduction, fast contact resolution, large deformation, St Venant tensor simulation, etc.) for future publications, although the PMI framework is equally applicable to them as well (e.g. Huang and Lee, 2013).

There are two ways to describe such articulated multiple rigid bodies: (1) maximal coordinates, which is the SE(3) rigid body transformation of each rigid body with the relative motion constraints implemented by Lagrange multipliers or through some other means (e.g. Hadap, 2006; Weinstein et al., 2006); and (2) generalized coordinates, which are obtained by reducing the degree of freedom (DOF) of the total system by applying their relative motion constraints (e.g. Lee et al., 2012; Spong et al., 2006). The maximal coordinate approach is typically preferred for simulating complex systems (e.g. systems composed of a large number of rigid bodies with possible collisions among them), whereas the generalized coordinate approach is more convenient or imperative for some mechanical systems (e.g. flexible beams with nodal coordinates, rigid links with elastic tendon routing, etc). The rendering framework proposed in this paper allows for both maximal and generalized coordinates as well as their combination (see Section 5.2 and Section 5.3).

More specifically, the maximal coordinate approach simulates the articulated multiple rigid bodies as a combination of their constituent single rigid body, whose dynamics in SE(3) can be written using the following Newton–Euler dynamics:

$$\begin{array}{c}\mathrm{mẍ}+b_xẋ+d{\psi }_x^T+{\lambda }_x=f+f_c\end{array}$$

(1)

$$\begin{array}{c}J\stackrel{.}{\omega }+\omega \times \mathrm{J\omega }+b_w\omega +d{\psi }_w^T+{\lambda }_w=\tau +{\tau }_c\end{array}$$

(2)

$$\begin{array}{c}Ṙ=\mathrm{RS}(\omega )\end{array}$$

(3)

where $x\in$ E(3) is the position with mass $m>0$ , $b_x\in {\Re }^{3\times 3}$ is a positive semi-definite damping matrix, and $f\in {\Re }^3$ is for embedding the exogenous force (e.g. virtual coupling), all expressed in the inertial frame $\{O\}$ . $R\in$ SO(3) is the rotation matrix of the body frame $\{B\}$ with respect to $\{O\}$ , $\omega \in$ so(3) is the angular velocity, $J\in {\Re }^{3\times 3}$ is the (constant) moment of inertia, $b_w\in {\Re }^{3\times 3}$ is a positive semi-definite damping matrix, and $\tau$ is for embedding the exogenous torque, all expressed in body frame $\{B\}$ . Here, ${\lambda }_x,{\lambda }_w\in {\Re }^3$ are to be used to enforce the constraints among the rigid bodies to realize their articulation, for which we utilize a potential-based approach (see Section 3.2); $f_c,{\tau }_c\in {\Re }^3$ are the contact forces such that $f_c^T\cdot ẋ+{\tau }_c^T\cdot w\le 0$ (e.g. with Coulomb friction; see Section 3.3); $d{\psi }_x,d{\psi }_w\in {\Re }^{1\times 3}$ are the one-form of certain potential energy $\psi (x,R)$ such that $\mathrm{d\psi }(x,R)∕\mathrm{dt}=d{\psi }_x\cdot ẋ+d{\psi }_w\cdot w$ (e.g. gravity or spring force; see Section 3.1); and $S(\star )$ is the skew-symmetric operator with $S(a)b=a\times b$ , $a,b\in {\Re }^3$ .

On the other hand, the dynamics of articulated rigid bodies in generalized coordinates can be written using the following continuous-time nonlinear Lagrangian dynamics:

$$\begin{array}{c}M(q)\stackrel{··}{q}+C(q,\stackrel{.}{q})\stackrel{.}{q}+B\stackrel{.}{q}+d{\psi }^T(q)=f+f_c\hfill \end{array}$$

(4)

where $q\in {\Re }^n$ is the configuration (i.e. generalized coordinates), $M(q)\in {\Re }^{n\times n}$ is the symmetric positive-definite inertia matrix, $C(q,\stackrel{.}{q})\in {\Re }^{n\times n}$ is the Coriolis matrix, $B\in {\Re }^{n\times n}$ is a positive semi-definite damping matrix, $\mathrm{d\psi }(q)\in {\Re }^{1\times n}$ is the one-form of potential energy function $\psi (q)\in \Re$ (e.g. gravity, spring potential) with $\mathrm{d\phi }∕\mathrm{dt}=\mathrm{d\phi }\cdot \stackrel{.}{q}$ , $f_c\in {\Re }^n$ is the contact force with $f_c^T\stackrel{.}{q}\le 0$ (e.g. with Coulomb friction), and $f\in {\Re }^n$ is for embedding the exogenous force (e.g. virtual coupling with a haptic device).

It is then well known that the mechanical system represented by maximal coordinates (1)–(3) or by generalized coordinates (4) satisfies the (continuous-time) passivity (Spong et al., 2006)

$$\begin{array}{c}{\int }_0^T{u}^{T}y\mathrm{dt}\ge E(T)-E(0),\forall T\ge 0\hfill \end{array}$$

(5)

with $E(t):=\kappa (t)+\psi (t)$ being the total energy, where $(u,y,\kappa )=([f;\tau ],[ẋ;w],(1∕2)m{ẋ}^Tẋ+(1∕2)w^T\mathrm{Jw})$ or $(u,y,\kappa )=(f,\stackrel{.}{q},(1∕2){\stackrel{.}{q}}^{T}M(q)\stackrel{.}{q})$ for the maximal coordinate formulation (1)–(3) or the generalized coordinate formulation (4), respectively. This passivity property (5) has been the key enabler of numerous nonlinear control techniques for mechanical systems, particularly of those exploiting their nonlinear dynamics (1)–(3) or (4) (i.e. passivity-based control; Lee, 2010; Lee and Lui, 2017; Slotine and Li, 1987; Spong et al., 2006) instead of unnecessarily attempting to fully eliminate it (e.g. feedback linearization).

Further, from the passivity theorem (Vidyasagar, 1993) that the feedback interconnection of any passive systems is necessarily stable regardless of their parameters (with some detectability condition met), this also implies that, if we can simulate each virtual entity to be passive, we may then construct very a complex virtual world simply by combining those passive virtual entities, each of which can be designed individually without referring to other components and whose parameters can be set arbitrarily as long as passivity is preserved (e.g. any positive-definite $M,B,K$ for (4), if $M(q)=M$ and $d{\psi }^T=\mathrm{Kq}$ ; see Section 2.2), while also theoretically ensuring the stability of the whole virtual world. The benefits of passive simulation, namely, modular design, parameter-invariant stability, and a theoretical stability guarantee, were recognized in Brown and Colgate (1997, 1998) and suggested as one of the ultimate goals of haptic rendering and interactive simulation.

This (fundamental) continuous-time passivity (5) is in general not yet easily achievable in the discrete-time domain, particularly when the simulation speed needs to be haptically fast. In fact, the work of Brown and Colgate (1998) shows that no explicit discrete-time passive integrators exist. To address this, in Lee and Huang (2008) we proposed a new integrator, NPMI, which is implicit, yet requires no iterations, and thus can still be solved haptically fast while also enforcing discrete-time passivity. This NPMI turns out to be a midpoint integrator, and is to be generalized in this paper to be a PMI. In Section 2.2, we summarize the result of Lee and Huang (2008), which stands as the backbone of the ensuing development of the PMI-based haptic rendering and interactive simulation framework of this paper.

2.2. Review of NPMI for linear mechanical systems

To introduce the main idea of NPMI (Lee and Huang, 2008), let us consider the following simple scalar linear mass–damper–spring dynamics:

$$\mathrm{ma}+\mathrm{bv}+\mathrm{kx}=f$$

where $x\in \Re$ , $v=ẋ\in \Re$ , and $a=ẍ\in \Re$ are the position, velocity and acceleration, $m,b,k>0$ are the mass, damper, and spring coefficients, and $f\in \Re$ is the exogenous force, respectively. The passivity property (5) is then given by the following. During $T_k=[t_k,t_{k+1})$ ,

$$\begin{array}{c}{\int }_{t_k}^{t_{k+1}}\mathrm{fv}\mathrm{dt}\hfill \\ =\frac{1}{2}m(v_{k+1}^2-v_k^2)+\frac{1}{2}k(x_{k+1}^2-x_k^2)+{\int }_{t_k}^{t_{k+1}}b{v}^2\mathrm{dt}\end{array}$$

where ${\star }_k:=\star (t_k)$ . Let us denote by ${\widehat{\star }}_k$ the “representative” ${\star }_k$ for the duration $T_k$ . Then, we can approximate the above passivity relation:

$$\begin{array}{c}{\widehat{f}}_k{\widehat{v}}_k{T}_k=\frac{1}{2}m{v}_{k+1}^2-\frac{1}{2}m{v}_k^2+\frac{1}{2}k{x}_{k+1}^2-\frac{1}{2}k{x}_k^2+b{\widehat{v}}_k^2{T}_k\hfill \\ =\frac{1}{2}(v_{k+1}-v_k)m(v_{k+1}+v_k)\\ +\frac{1}{2}(x_{k+1}-x_k)k(x_{k+1}+x_k)+b{\widehat{v}}_k^2{T}_k\end{array}$$

(6)

where we also use $T_k:=t_{k+1}-t_k$ with a slight abuse of notation.

Our goal here is then to find a discrete-time integrator map

$$\mathcal{L}:(x_k,v_k,f_k)\to (x_{k+1},v_{k+1})$$

and the representative power variables ${\widehat{v}}_k,{\widehat{f}}_k$ so that the map $\mathcal{L}$ exactly duplicates the continuous-time passivity relation (6), or, equivalently, its differential form

$$\begin{array}{c}{f}^{T}v=v(\mathrm{ma}+\mathrm{bv}+\mathrm{kx})\hfill \end{array}$$

(7)

By comparing (6) with (7), the NPMI then suggests the following kinematics relation (Lee and Huang, 2008):

$$\begin{array}{c}{\widehat{v}}_k:=\frac{v_{k+1}+v_k}{2}=\frac{x_{k+1}-x_k}{T_k}\hfill \end{array}$$

(8)

and, further, applying (8) to (6) and comparing with (7), it also suggests the following dynamics update equation:

$$\begin{array}{c}m\frac{v_{k+1}-v_k}{T_k}+b{\widehat{v}}_k+k{\widehat{x}}_k=f_k\hfill \end{array}$$

(9)

with ${\widehat{x}}_k:=\frac{x_{k+1}+x_k}{2}$ and ${\widehat{f}}_k=f_k$ .

This NPMI (8)–(9) then exactly duplicates the continuous-time passivity property (5) by using the “representative” velocity ${\widehat{v}}_k$ and position ${\widehat{x}}_k$ . For the choice of ${\widehat{f}}_k$ , we have latitude, and thus we choose ${\widehat{f}}_k=f_k$ for simplicity of the integrator $\mathcal{L}$ . With these “representative” quantities, the NPMI then satisfies the following discrete-time passivity:

$$\begin{array}{c}\sum _{k=0}^{N-1}{f}_k{\widehat{v}}_k{T}_k=\frac{1}{2}m(v_N^2-v_o^2)+\frac{1}{2}k(x_N^2-x_o^2)+\sum _{k=0}^{N-1}b{\widehat{v}}_k^2{T}_k\hfill \\ \ge -\frac{1}{2}m{v}_o^2-\frac{1}{2}k{x}_o^2=-E_o\end{array}$$

(10)

$\forall N\ge 0$ , in other words, the maximum energy extractable from the system is lower-bounded by the system’s initial energy. Since it is derived directly by enforcing a certain physical principle (i.e. passivity) in the discrete-time domain, the NPMI is a geometric integrator similar to variational integrators (e.g. Marsden and West, 2001), which discretizes the variational principle itself.

This NPMI (8)–(9) also turns out to be a midpoint integrator (e.g. Bathe, 2007; Krysl, 2005, 2006; Kuhl and Crisfield, 1999; Simo et al., 1992), which is well known for its superior energy-conserving property (i.e. passivity, but with no input–output consideration) in the structural mechanics community, but has not been adopted so far for haptic rendering and interactive simulation. In this paper, we generalize this NPMI as a PMI, with full consideration of input–output mapping (i.e. power port), generalized and maximal coordinate formalisms, nonlinear dynamics (1)–(4), and constraints and multi-point contacts with Coulomb friction, all the while retaining its discrete-time passivity and non-iterativeness by utilizing certain coordinate transformations and potential action approximation.

Of course, other types of integrators have also been used for haptic rendering and interactive simulation, namely, explicit Euler integrators (EEIs), symplectic Euler integrators (SEIs), and implicit Euler integrators (IEIs), with the second being adopted for the majority of widely used open physics engines (e.g. Bullet and ODE, with VREP and CHAI-3D built upon them). In Section 2.3, we perform a simple yet effective comparison study among EEI, SEI, IEI and PMI.

2.3. Comparison of EEI, SEI, IEI and PMI

Consider the simple scalar linear mass–spring system, which can be written as

$$\begin{array}{cc}m\frac{v_{k+1}-v_k}{T_k}+k{x}_k=f_k\hfill & x_{k+1}=x_k+v_k{T}_k\\ m\frac{v_{k+1}-v_k}{T_k}+k{x}_k=f_k\hfill & x_{k+1}=x_k+v_{k+1}{T}_k\\ m\frac{v_{k+1}-v_k}{T_k}+k{x}_{k+1}=f_k\hfill & x_{k+1}=x_k+v_{k+1}{T}_k\\ m\frac{v_{k+1}-v_k}{T_k}+k{\widehat{x}}_k=f_k\hfill & x_{k+1}=x_k+{\widehat{v}}_k{T}_k\end{array}$$

presented in the order EEI, SEI, IEI and PMI, where ${\widehat{x}}_k,{\widehat{v}}_k$ are defined in (8) and (9). Note the difference in their rendering spring force: $k{x}_k$ (EEI and SEI), $k{x}_{k+1}$ (IEI), and $k{\widehat{x}}_k$ (PMI). Note also that each integrator adopts a different “representative” velocity: $v_k$ (EEI), $v_{k+1}$ (SEI and IEI), and ${\widehat{v}}_k$ (PMI).

Then, adopting a procedure similar to that in Section 2.2, we can obtain the passivity property of each integrator as follows. First, for the EEI, we have

$$\sum _{k=0}^{N-1}{f}_k{v}_k{T}_k\le \frac{1}{2}m(v_N^2-v_0^2)+\frac{1}{2}k(x_N^2-x_0^2)$$

for which we utilize

$$\begin{array}{c}\sum _{k=1}^{N}m\frac{v_{k+1}-v_k}{T_k}{v}_k{T}_k\le \frac{1}{2}m{v}_{N+1}^2-\frac{1}{2}m{v}_0^2\\ \sum _{k=1}^{N}k{x}_k\frac{x_{k+1}-x_k}{T_k}{T}_k\le \frac{1}{2}k{x}_{N+1}^2-\frac{1}{2}k{x}_0^2\end{array}$$

(11)

derived by using the Cauchy–Schwarz inequality. This then implies that EEI is in general not passive. Similarly, for the IEI, we have

$$\sum _{k=0}^{N-1}{f}_k{v}_{k+1}{T}_k\ge \frac{1}{2}m(v_N^2-v_0^2)+\frac{1}{2}k(x_N^2-x_0^2)$$

for which we utilize

$$\begin{array}{c}\sum _{k=1}^{N}m\frac{v_{k+1}-v_k}{T_k}{v}_{k+1}{T}_k\ge \frac{1}{2}m{v}_{N+1}^2-\frac{1}{2}m{v}_0^2\\ \sum _{k=1}^{N}k{x}_{k+1}\frac{x_{k+1}-x_k}{T_k}{T}_k\ge \frac{1}{2}k{x}_{N+1}^2-\frac{1}{2}k{x}_0^2\end{array}$$

(12)

This then implies that IEI is passive, but in general dissipative, and thus not able to simulate lossless behavior such as harmonic oscillation. On the other hand, for SEI, its passivity is undetermined, since its relation is a combination of (12) (i.e. passive) and (11) (i.e. non-passive). This in fact matches well with the known fact that SEI is only conditionally stable depending on parameters. Finally, for PMI, we have

$$\sum _{k=0}^{N-1}{f}_k{\widehat{v}}_k{T}_k=\frac{1}{2}m(v_N^2-v_0^2)+\frac{1}{2}k(x_N^2-x_0^2)$$

from (10), implying that it is passive, and, furthermore, exactly lossless, and thus can precisely simulate harmonic oscillation.

This observation is then confirmed by the simulation results in Figure 1 and Figure 2, where (1) the EEI always diverges; (2) the SEI can simulate harmonic oscillation with large mass, yet diverges with small mass; (3) the IEI always converges to zero, yet is not able to sustain harmonic oscillation, clearly showing its dissipativeness; and (4) only the PMI can simulate harmonic oscillation with the total energy exactly conserved for any mass/spring ratio. This then clearly shows the advantage of our PMI, namely its usability for passive (i.e. stable) simulation with a wide range of parameters as compared to EEI and SEI, while also being typically simpler to implement than IEI on top of its ability to simulate lossless behavior.

Fig. 1.

Simple harmonic oscillation simulation with EEI, SEI, IEI and PMI: $m=0.01$ kg, $k=100$ N/m, $T_k=0.001$ s.

Fig. 2.

Simple harmonic oscillation simulation with EEI, SEI, IEI and PMI: $m=0.00001$ kg, $k=100$ N/m, $T_k=0.001$ s.

3. PMI-based simulation in maximal coordinates

The NPMI (8)–(9) of Section 2.2 exactly enforces passivity only for linear mechanical systems. Here, we generalize this NPMI to PMI of rigid bodies in SE(3), which constitute the basic building block of the maximal coordinate formalism of articulated rigid bodies. We also show how to incorporate constraints into this PMI-based simulation of rigid bodies to realize articulation among the rigid bodies, and how to include multi-point Coulomb friction contact via LCP formulation into the framework of PMI. The main complication with this is that the configuration space of rigid bodies in maximal coordinates is SE(3), which is not a vector space. The multi-point contact computation problem is also typically complex and time-consuming. To circumvent these challenges, here we utilize some approximations in enforcing passivity. Such approximation is in general necessary for other haptic rendering or interactive simulation results to attain interactivity (e.g. contact force linearization with IEI or SEI; Otaduy and Lin, 2006) or for structural dynamics to speed up the computation (e.g. interaction explicitation; Bathe, 2007; Krysl, 2005, 2006), which is known to work adequately in most cases, as confirmed through our experiments as well (see Section 5).

3.1. PMI of rigid body in SE(3)

Consider the Newton–Euler dynamics (1)–(3) of rigid bodies in SE(3). In this section, we derive the PMI formulation of (1)–(3). For this, we consider only the inertial dynamics (i.e. associated with $m,J$ ), the potential action (i.e. $d{\psi }_{\star }$ ), and the exogenous forcing $f,\tau$ . The constraint and contact forces (i.e. ${\lambda }_{\star },f_c,{\tau }_c$ ) will be added in Section 3.2 and Section 3.3, respectively. We also restrict the potential $\psi$ in this section to be the linear E(3) spring directly acting on x (i.e. $\psi ={\psi }_x=\frac{1}{2}{x}^T{k}_{x}x$ ) or the SO(3) spring acting on R (i.e. $\psi ={\psi }_w(R)$ ), while postponing the treatment of more general $\psi (x,R)$ to Section 3.2.

Then, with the linear spring potential ${\psi }_x=\frac{1}{2}{x}^T{k}_{x}x$ , the translation dynamics (1) can simply be simulated via the linear NPMI (8)–(9) while exactly enforcing passivity as stated in Section 2.2. The attitude dynamics (2) also has a structure similar to that of linear NPMI (9) except for the coupling term $w\times \mathrm{Jw}$ , which can be written as $w\times \mathrm{Jw}=-S(\mathrm{Jw})w$ , where $S(\mathrm{Jw})\in {\Re }^{3\times 3}$ is skew-symmetric with $a^{T}S(\mathrm{Jw})a=0\forall a\in {\Re }^3$ . This then suggests the following PMI expression of the attitude dynamics (2):

$$\begin{array}{c}J\frac{{\omega }_{k+1}-{\omega }_k}{T_k}-S(J{\omega }_k){\widehat{\omega }}_k+b_w{\widehat{\omega }}_k+d{\psi }_{w_k}^T={\tau }_k\hfill \end{array}$$

(13)

with the “representative” angular velocity

$$\begin{array}{c}{ŵ}_k:=\frac{w_{k+1}+w_k}{2}\hfill \end{array}$$

(14)

similar to (8). It is then not difficult to show the following passivity (lossless) property with ${\psi }_w=0$ similar to (10):

$$\begin{array}{c}\sum _{k=0}^{N-1}{\tau }_k^T{\widehat{\omega }}_k{T}_k=\frac{1}{2}{w}_N^{T}J{w}_N-\frac{1}{2}{w}_o^{T}J{w}_o\ge -\frac{1}{2}{w}_o^{T}J{w}_o\hfill \end{array}$$

(15)

where we use ${ŵ}_k^{T}S(J{w}_k){ŵ}_k=0$ , $\forall k=0,\dots ,N-1$ . The coupling term $S(\mathrm{Jw})w$ is often omitted in computer graphics literature. For this case, the PMI expression (13) simply reduces to the linear NPMI (9).

One thing still missing here is how to connect the angular velocity ${ŵ}_k\in {\Re }^3$ in (14) to the rotation matrix $R_k\in$ SO(3) while ensuring passivity as much as is permissible for the interactivity requirement. For this, we utilize the following discrete-time version of the rotational kinematics (3):

$$\begin{array}{c}{R}_{k+1}=R_{k}exp(S({\widehat{\omega }}_k{T}_k))\hfill \end{array}$$

(16)

where $\exp (\star )$ is the exponential map (Murray et al., 1993). This discrete-time kinematics equation (16) in SO(3) then should be compatible with the SO(3) spring potential, as so are true for ${\widehat{v}}_k,{\widehat{x}}_k$ and ${\psi }_k=\frac{1}{2}k{x}_k^2$ in Section 2.2. For this, define the following SO(3) spring potential:

$$\begin{array}{c}{\psi }_w:=\frac{1}{2}{k}_w{\varphi }^2(R,R_o)\hfill \end{array}$$

(17)

with $k_w>0$ and

$$\begin{array}{c}\varphi (R,R_o):=\parallel \log (R^T{R}_o)\parallel =:{\cos }^{-1}\gamma \hfill \end{array}$$

where $R_o\in$ SO(3) is the reference rotation, $\mathrm{tr}(\star )$ and $\log (\star )$ are the trace and matrix logarithm operators, and

$$\gamma =\frac{1}{2}[\mathrm{tr}(R^T{R}_o)-1]$$

This $\varphi$ is known to be the geodesic on the unit sphere (Huynh, 2009). This SO(3) spring, of course, is practically important as well, for example directly applicable to implementing virtual coupling between two large-rotation rigid bodies (e.g. virtual coupling for palm and fingers of a multi-fingered hand; see Figure 4).

Our objective is then to find the action $d{\psi }_{w_k}$ for (13) to implement the SO(3) spring (17) under the passivity requirement. For this, particularly since SO(3) is not a vector space unlike for the case of the linear NPMI (9), we attempt to enforce passivity only approximately. This, however, turns out to work fairly well as shown below. More precisely, we aim to define $d{\psi }_{w_k}$ such that

$$\begin{array}{c}d{\psi }_{w_k}{\widehat{\omega }}_k{T}_k=\frac{\partial {\psi }_w}{\mathrm{\partial \varphi }}\Delta {\varphi }_k+\frac{1}{2}\frac{{\partial }^2{\psi }_w}{\partial {\varphi }^2}\Delta {\varphi }_k^2\hfill \end{array}$$

(18)

which is the second-order approximation of the exact potential difference (i.e. ${\psi }_{w,k+1}-{\psi }_{w_k}$ ), with $\partial {\psi }_w∕\mathrm{\partial \varphi }=k_w{\varphi }_k$ and ${\partial }^2{\psi }_w∕\partial {\varphi }^2=k_w$ . For (18), we also use

$$\begin{array}{c}\Delta {\varphi }_k=\frac{\mathrm{\partial \varphi }}{\mathrm{\partial \gamma }}\Delta {\gamma }_k,\Delta {\gamma }_k:=\frac{1}{2}{\widehat{\omega }}_k^T{\bar{\omega }}_k{T}_k\hfill \end{array}$$

(19)

where $\mathrm{\partial \varphi }∕\mathrm{\partial \gamma }=-1∕\sin {\varphi }_k$ and ${\bar{\omega }}_k\in {\Re }^3$ is defined such that

$$S({\bar{\omega }}_k):=R_k^T{R}_o-R_o^T{R}_k$$

which has the direction of the equivalent axis between $R_k$ and $R_o$ and the magnitude $\parallel {\bar{\omega }}_k\parallel =2\sin ({\varphi }_k)$ . The term $\Delta {\gamma }_k$ is designed such that

$$\begin{array}{c}\Delta {\gamma }_k=\frac{1}{2}\mathrm{tr}[R_{k+1}^T{R}_o-R_k^T{R}_o]\approx \frac{1}{2}\mathrm{tr}[S^T({\widehat{\omega }}_k{T}_k)R_k^T{R}_o]\hfill \\ =\frac{1}{4}\mathrm{tr}[S^T({\widehat{\omega }}_k{T}_k)(R_k^T{R}_o-R_o^T{R}_k)]\end{array}$$

which is reduced to (19), where we use (16), the first-order approximation of $\exp (S({ŵ}_k{T}_k))\approx I+S({ŵ}_k{T}_k)$ , and the properties of $\mathrm{tr}(\star )$ , for example, $\mathrm{tr}(A+B)=\mathrm{tr}(A)+\mathrm{tr}(B)$ , $\mathrm{tr}(\mathrm{ABC})=\mathrm{tr}(\mathrm{BCA})$ , $\mathrm{tr}(A)=\mathrm{tr}(A^T)$ , and $a^{T}b=\mathrm{tr}[S^T(a)S(b)∕2]$ , $\forall a,b\in {\Re }^3$ .

Applying these to (18), we can then define $d{\psi }_{w_k}$ :

$$\begin{array}{c}d{\psi }_{w_k}^T=-\frac{k_w{\varphi }_k}{2\sin {\varphi }_k}{\bar{\omega }}_k+\frac{1}{8}\frac{k_w{T}_k}{{\sin }^2{\varphi }_k}({\bar{\omega }}_k{\bar{\omega }}_k^T){\widehat{\omega }}_k\hfill \end{array}$$

where ${\bar{\omega }}_k{\bar{\omega }}_k^T\in {\Re }^{3\times 3}$ is a rank-one matrix. Since $\left|\right|{\bar{\omega }}_k\left|\right|=2\sin ({\varphi }_k)$ , we can further approximately simplify $d{\psi }_{w_k}$ :

$$\begin{array}{c}d{\psi }_{w_k}^T\approx -\frac{k_w{\varphi }_k}{2\sin {\varphi }_k}{\bar{\omega }}_k+\frac{k_w{T}_k}{2}{\widehat{\omega }}_k\hfill \end{array}$$

(20)

Note that the second term acts similarly to viscous damping action. The passivity of this SO(3) spring action (20) is then demonstrated in Figure 3, where stable sustained oscillation is attained with the total energy almost precisely preserved even if the action (20) is designed to enforce passivity only approximately.

Fig. 3.

Total energy of rotational dynamics (13) with SO(3) spring (20) with nearly perfect energy conservation and stable/sustained oscillation in SO(3): with $J=0.001\cdot I$ kg m², $k_w=1$ Nm/rad and $T_k=0.01$ s starting from ${\varphi }_o=0.9\pi$ along the direction of $\frac{1}{\sqrt{6}}(1,2,1)$ .

Finally, the PMI expression of the Newton–Euler dynamics (1)–(2) with the linear spring for x and the SO(3) spring potential (20) can be written in one single equation:

$$M\frac{V_{k+1}-V_k}{T_k}+C_k({\omega }_k){\widehat{V}}_k+B{\widehat{V}}_k+d{\psi }_k^T=F_k$$

where $V_k:=[v_k;{\omega }_k]\in {\Re }^6$ with ${\widehat{V}}_k:=[{\widehat{v}}_k;{\widehat{\omega }}_k]\in {\Re }^6$ , $M:=\mathrm{tr}[m{I}_{3\times 3},J]\in {\Re }^{6\times 6}$ is the inertia matrix, $C_k:=\mathrm{tr}[0_{3\times 3},-S(J{\omega }_k)]$ with $C_k=-C_k^T$ , $B:=\mathrm{tr}[b_x,b_{\omega }]\in {\Re }^{6\times 6}$ is the damping matrix, $d{\psi }_k^T:=[k_x{\widehat{x}}_k;d{\psi }_{{\omega }_k}^T]\in {\Re }^6$ is the potential action as defined in (9) and (20), and $F_k:=[f_k;{\tau }_k]\in {\Re }^6$ is the exogenous wrench. This equation can also be written in the following form, which is to be used in the ensuing developments:

$$\begin{array}{c}{\widehat{M}}_k{\widehat{V}}_k=M_k{V}_k-d{\bar{\psi }}_k^T+F_k\hfill \end{array}$$

(21)

where ${\widehat{M}}_k:=[M_k+C_k+B+(T_k∕2)K]\in {\Re }^{6\times 6}$ with $M_k:=2M∕T_k$ , $K:=\mathrm{tr}[k_x,k_{\omega }]\in {\Re }^{6\times 6}$ is the spring matrix from (9) and (20), and $d{\bar{\psi }}_k^T:=[-k_x{x}_k;\frac{k_{\omega }{\varphi }_k}{2\sin {\varphi }_k}{\bar{\omega }}_k]\in {\Re }^6$ from (20). Here, both ${\widehat{M}}_k$ and $M_k$ are non-singular.

3.2. PMI-based rendering of articulated rigid bodies

Some of the most widely used entities in haptic rendering and interactive simulation are so-called articulated rigid bodies (e.g. Weinstein et al., 2006). Here, we propose to construct such articulated rigid bodies as a collection of multiple rigid links with their joints implemented by constraining the positions of some points on successive links via very stiff springs. For instance, see the multi-fingered hand in Figure 4, where the distal interphalangeal (DIP) and metacarpophalangeal (MCP) joints of the index finger are implemented by enforcing two-point position constraints on the distal-middle and middle-proximal phalanges, whereas the thumb carpometacarpal (CMC) joint is implemented by one single-point position constraint between the palm and the metacarpal bone. Other types of joint constraints can also be implemented in a similar manner by using different spring potentials (e.g. elliptic joint with four-point springs (Knoop et al., 2005); prismatic joint with an SO(3) potential and a two-directional translation spring potential, etc.), details of which will be reported in future publications.

Fig. 4.

Multi-fingered hand rendered as articulated multiple rigid bodies with joints implemented by point constraint springs.

Fig. 5.

Snapshots of rigid box falling on rigid level surface: it bounces back and forth under PMI-LCP (white box), but completely plastically collides under SEI (black box) even with the same parameter settings.

Consider m rigid links with l point position constraints (e.g. m=10 and l=16 for the hand in Figure 4). We assume each of these l points is shared only by two successive links, although the development presented here can easily be extended to more general cases. Then, we can define the following map:

$$c_1,c_2:\{1,2,\dots ,l\}\mapsto \{1,2,\dots ,m\}$$

in other words, connecting the i th position constraining point ( $i\in \{1,2,\dots ,l\}$ ) to its presence on the $c_1(i)$ th and the $c_2(i)$ th links. We can then denote this presence of the i th constraining point as

$$p_{i,j}:=x_{c_j(i)}+R_{c_j(i)}{r}_{c_j(i),i}\in {\Re }^3,j=1,2$$

which is on the $c_j(i)$ th link with the offset from $x_{c_j(i)}$ given by $r_{c_j(i),i}\in {\Re }^3$ expressed in the $c_j(i)$ th link body fixed frame.

To implement the joints among the m links, we then want to enforce

$$\Delta p_i:=p_{i,1}-p_{i,2}\equiv 0$$

Stacking up all the l point position constraints, define the constraint function $G:{\Re }^{3l}\to {\Re }^{3l}$ as

$$G:=[\Delta p_1;\Delta p_2;\dots \Delta p_l]$$

with $G=0$ implying that all the l point position constraints are ensured, and so are the joints among the m links. We can then first-order approximate this G:

$$\begin{array}{c}{G}_{k+1}\approx G_k+J_k{\widehat{V}}_k{T}_k\hfill \end{array}$$

(22)

where ${\widehat{V}}_k=[{\widehat{v}}_{1_k};{\widehat{\omega }}_{1_k};{\widehat{v}}_{2_k};{\widehat{\omega }}_{2_k};\dots ;{\widehat{v}}_{m_k};{\widehat{\omega }}_{m_k}]\in {\Re }^{6m}$ contains the representative velocities of all the m rigid bodies, and $J\in {\Re }^{3l\times 6m}$ is the Jacobian defined such that $\mathrm{dG}∕\mathrm{dt}=\mathrm{JV}$ .

The PMI expression of the dynamics of the m rigid links with the l constraints can then be written similarly to (21):

$$\begin{array}{c}{\widehat{M}}_k\widehat{V_k}=M_k{V}_k-d{\bar{\psi }}_k^T+J_k^T{\lambda }_k+F_k\hfill \end{array}$$

(23)

where ${\lambda }_k\in {\Re }^{3l}$ is for embedding certain potential action to robustly maintain the constraint G=0, similar to the penalty-based contact formulation (Drumwright, 2008; Hasegawa and Fujii, 2003), and all the other terms are defined for the m rigid bodies similarly to those in (21). Our goal is then to design the potential action ${\lambda }_k$ to attain the constraint $G=0$ while also respecting the passivity requirement. For this, let us define the following potential function associated with the (target) constraint G=0 :

$$\begin{array}{c}{\psi }_{{\lambda }_k}:=\frac{1}{2}{k}_{\lambda }{G}_k^T{G}_k\hfill \end{array}$$

where $k_{\lambda }>0$ is a constant.

Then, using (10), (15), (18), (22) and (23), we can obtain

$$\begin{array}{c}{F}_k^T{\widehat{V}}_k{T}_k\approx {\kappa }_{k+1}-{\kappa }_k+{\psi }_{k+1}-{\psi }_k+{\widehat{V}}_k^{T}B{\widehat{V}}_k{T}_k-{\widehat{V}}_k^T{J}_k^T{\lambda }_k{T}_k\end{array}$$

where ${\kappa }_k:=\frac{1}{2}{V}_k^{T}M{V}_k$ is the kinetic energy, M is the inertia matrix and ${\psi }_k=\frac{1}{2}{x}_k^T{k}_x{x}_k+{\psi }_{w,k}$ is the potential energy of the m rigid bodies collective as defined before (21). Further, utilizing the second-order approximation of

$$\begin{array}{c}{\psi }_{{\lambda }_{k+1}}-{\psi }_{{\lambda }_k}\approx k_{\lambda }{G}_k^T{J}_k{\widehat{V}}_k{T}_k+\frac{1}{2}{k}_{\lambda }{\widehat{V}}_k^T{J}_k^T{J}_k{\widehat{V}}_k{T}_k^2\hfill \end{array}$$

(24)

we can design

$$\begin{array}{c}{\lambda }_k=-k_{\lambda }[G_k+\frac{1}{2}{J}_k{\widehat{V}}_k{T}_k]\hfill \end{array}$$

(25)

with which the above relation reduces to the following: with B=0,

$$\begin{array}{cc}{F}_k^T{\widehat{V}}_k{T}_k\hfill & \approx E_{k+1}-E_k\end{array}$$

where $E_k:=\frac{1}{2}{V}_k^{T}M{V}_k+{\psi }_k+{\psi }_{\lambda ,k}$ is the total energy. With this ${\lambda }_k$ from (25), (23) can also be written as

$$\begin{array}{cc}\left[{\widehat{M}}_k+\frac{k_{\lambda }{T}_k}{2}{J}_k^T{J}_k\right]\hfill & {\widehat{V}}_k=M_k{V}_k-k_{\lambda }{J}_k^T{G}_k-d{\bar{\psi }}_k^T+F_k\end{array}$$

(26)

The analysis above shows that the potential action ${\lambda }_k$ (25) would robustly enforce the constraint G=0 while also nearly enforcing the passivity of the PMI simulation of articulated rigid bodies. Here, note that the passivity is enforced only approximately, as we adopt the potential action approximation (24), which in turn stems from our adoption of the Jacobian approximation (22). We have this nonlinear potential approximation (24) even if the potential ${\psi }_{{\lambda }_k}$ itself is only quadratic, since the maximal coordinates are on a non-vector space SE(3). This potential action approximation (24) is crucial to rendering the PMI simulation still non-iterative even with the nonlinear potential ${\psi }_{{\lambda }_k}$ . The “lack” of passivity due to this potential action approximation (24) can also be easily amended by injecting a (typically) small amount of damping B, which can also be incorporated into the PMI simulation while retaining its non-iterative structure (26).

This passive enforcing of the constraint $G=0$ turns out to be very useful, since it allows us to utilize very large stiffness (i.e. $k_{\lambda }$ ) to enforce the joint articulations, even for very light rigid links (e.g. to render light the virtual proxy for transparency); see Section 5.1. Note that we can also implement the spring potential action $\psi (x,R)$ (e.g. spring connected not on the center of mass x of rigid bodies) by using the same procedure to design ${\lambda }_k$ proposed here, that is, define the spring connection point positions similarly to G, define their spring potentials similarly to ${\psi }_{\lambda }$ , and design the potential action $\mathrm{d\psi }$ similarly to ${\lambda }_k$ to respect the approximated potential difference equation. Here, we also adopt this “soft” potential action ${\psi }_{{\lambda }_k}$ to enforce the constraint G=0 instead of the “hard” Lagrange multiplier approach for (22), since (1) the latter is much less robust (e.g. drift) with the Jacobian approximation (22) than the former, which contains potential action similar to the penalty-based contact formulation (Drumwright, 2008; Hasegawa and Fujii, 2003); and (2) the passivity of the former is easily achievable (with small damping B injection), whereas that of the latter still unclear.

3.3. PMI-LCP for multi-point Coulomb friction contact

The contact problem has been one of the key problems in computer graphics, structural mechanics, and haptic rendering. For this, there are two main directions of research, namely the constraint-based method (e.g. Duriez et al., 2006; Peterlik et al., 2011) and the penalty-based method (e.g. Drumwright, 2008; Kolesnikov and Zefran, 2007). The penalty-based method generates repulsive elastic force when two objects overlap with each other. This method is simple and easy to implement, but it has inherent limitations: allowing violation of the contact constraint (e.g. penetration). It is in general not straightforward either to implement Coulomb friction force with this penalty-based method. On the other hand, the constraint-based method directly calculates contact force, which is necessary for satisfying the contact constraint. By doing so, this method prevents violation of the constraint (to the extent of collision detection error and numerical precision error though, which can induce constraint drift). Friction generation is also more convenient with the constraint-based method.

To render the contact among rigid objects, we adopt the point contact model with Coulomb friction. Consider a contact point between two colliding rigid bodies. We can then define the contact normal direction $\overrightarrow{n}(t)\in {\Re }^3$ and the penetration depth along that direction. Denote this depth by $d(t)\ge 0$ , with $d(t)=0$ and $d(t)>0$ respectively implying contact and no contact. Whenever $d(t)=0$ occurs, the following higher-order conditions can be utilized to enforce the no-penetration constraint:

$$\begin{array}{c}ḋ\ge 0,\stackrel{··}{d}\ge 0,\dots \hfill \end{array}$$

with the non-negative signs allowing only for separation or maintaining the contact. In this paper, we use the velocity-level constraint, which can then be written as the following complementarity condition (Siciliano and Khatib, 2008):

$$\begin{array}{c}0\le ḋ\perp {\lambda }^n\ge 0\hfill \end{array}$$

(27)

where ${\lambda }^n\ge 0$ is the contact force along $\overrightarrow{n}$ .

Now, suppose that the articulated rigid bodies (23) interact with another virtual rigid object, and denote the i th contact point from that contact by

$$p_i^c:=x_{c(i)}+r_i^c$$

which is on the $c(i)$ th link and occurs at the offset of $r_c^i\in {\Re }^3$ from its center of mass (expressed in the inertial frame), where $c:\{1,2,\dots ,l_c\}\to \{1,2,\dots ,m\}$ is the mapping from the contact points to their associated rigid links from the articulated rigid bodies (23). The contact wrench, which is to be embedded in $F_k$ of (23), is then given by

$$F_{i,k}^c={\underset{=:A{d}_{i,k}^c\in {\Re }^{3\times 6}}{\underbrace{[\begin{array}{cc}I& -S(r_{i,k}^c)R_{\sigma (\mathrm{ }i),k}\end{array}]}}}^T{f}_{i,k}^c\in {\Re }^6$$

(28)

where $A{d}_{i,k}^c$ is the (partial) contact adjoint operator (Murray et al., 1993). The contact force $f_{i,k}^c\in {\Re }^3$ is defined following the point contact and Coulomb friction model (Siciliano and Khatib, 2008) such that

$$\begin{array}{c}{f}_{i,k}^c={\overrightarrow{n}}_k{\lambda }_{i,k}^n+\left[\begin{array}{c}{\overrightarrow{t}}_{i,k}^1,{\overrightarrow{t}}_{i,k}^2,\dots ,{\overrightarrow{t}}_{i,k}^{\bar{p}}\hfill \end{array}\right]\left[\begin{array}{c}{\lambda }_{i,k}^{t_1};{\lambda }_{i,k}^{t_2};\dots ;{\lambda }_{i,k}^{t_{\bar{p}}}\hfill \end{array}\right]\hfill \\ =:N_{i,k}{\lambda }_{i,k}^n+D_{i,k}{\lambda }_{i,k}^t\end{array}$$

(29)

where ${\lambda }_{i,k}^n\ge 0$ is the normal contact force, ${\lambda }_{i,k}^t\in {\Re }^{\bar{p}}$ are the tangential friction forces, for which we utilize the $\bar{p}$ -sides pyramidal cone approximation (Duriez et al., 2006), ${\overrightarrow{n}}_{i,k},{\overrightarrow{t}}_{i,k}^j\in {\Re }^3$ are the normal and tangential direction vectors, and $N_{i,k}\in {\Re }^3,D_{i,k}\in {\Re }^{3\times \bar{p}}$ are the stacking matrices of ${\overrightarrow{n}}_{i,k}$ and ${\overrightarrow{t}}_{i,k}^j$ .

The PMI expression of the articulated rigid bodies, the virtual object, and their contact interaction can then be obtained by augmenting (26):

$$\begin{array}{c}\underset{=:A_k\in {\Re }^{6(m+1)\times 6(m+1)}}{\underbrace{\left[\begin{array}{cc}{\widehat{M}}_k+\frac{k_{\lambda }{T}_k}{2}{J}_k^T{J}_k\hfill & 0\hfill \\ 0\hfill & {\widehat{M}}_k^o\hfill \end{array}\right]}}\underset{=:{\widehat{V}}_k^a\in {\Re }^{6(m+1)}}{\underbrace{\left(\begin{array}{c}{\widehat{V}}_k\hfill \\ {\widehat{V}}_k^o\hfill \end{array}\right)}}\\ =\underset{=:B_k\in {\Re }^{6(m+1)}}{\underbrace{\left(\begin{array}{c}{M}_k{V}_k-k_{\lambda }{J}_k^T{G}_k-d{\bar{\psi }}_k^T\hfill \\ M_k^o{V}_k^o\hfill \end{array}\right)}}+F_k^c\end{array}$$

with the contact force term $F_k^c\in {\Re }^{6(m+1)}$ defined by

$$F_k^c:={\bar{N}}_k{\lambda }_k^n+{\bar{D}}_k{\lambda }_k^t$$

where ${\lambda }_k^n:=[{\lambda }_{1,k}^n;{\lambda }_{2,k}^n;\dots ;{\lambda }_{l_c,k}^n]\in {\Re }^{l_c}$ and ${\lambda }_k^t:=[{\lambda }_{1,k}^t;{\lambda }_{2,k}^t;\dots ;{\lambda }_{l_c,k}^t]\in {\Re }^{l_c\times \bar{p}}$ are the collections of the normal contact and tangential friction forces of all $l_c$ contact points, and the matrices ${\bar{N}}_k\in {\Re }^{6(m+1)\times l_c}$ and ${\bar{D}}_k\in {\Re }^{6(m+1)\times (l_c\times \bar{p})}$ are the matrices specifying the directions of ${\lambda }_k^n$ and ${\lambda }_k^t$ with $A{d}_{i,k}^c$ in (28) embedded in them. We can then rewrite the above PMI simulation equation:

$$\begin{array}{c}{\widehat{V}}_k^a=A_k^{-1}(B_k+{\bar{N}}_k{\lambda }_k^n+{\bar{D}}_k{\lambda }_k^t)\hfill \end{array}$$

(30)

Similar to (27), here we adopt the velocity-level constraint. We can then write the following complementarity conditions (Lloyd, 2005):

$$\begin{array}{c}0\le {\sigma }_k^n\perp {\lambda }_k^n\ge 0\end{array}$$

(31)

$$\begin{array}{c}0\le {\sigma }_k^t\perp {\lambda }_k^t\ge 0\end{array}$$

(32)

$$\begin{array}{c}0\le {\sigma }_k^f\perp {\lambda }_k^f\ge 0\end{array}$$

(33)

where

$$\begin{array}{c}{\sigma }_k^n:={\bar{N}}_k^T{\widehat{V}}_k^a\in {\Re }^{l_c}\\ {\sigma }_k^t:={\bar{D}}_k^T{\widehat{V}}_k^a+U{\lambda }_k^f\in {\Re }^{l_c\times \bar{p}}\\ {\sigma }_k^f:=\bar{\mu }{\lambda }_k^n-U^T{\lambda }_k^t\in {\Re }^{l_c}\end{array}$$

with $\bar{\mu }:=\mu I_{l_c\times l_c}$ and $U:=\mathrm{tr}[1_{\bar{p}},1_{\bar{p}},\dots 1_{\bar{p}}]\in {\Re }^{(l_c\times \bar{p})\times \bar{p}}$ . Here, (31) and (32) respectively specify the complementarity conditions along the normal and tangential directions, whereas (33) the Coulomb friction model where the friction cannot exceed $\mu {\lambda }^n$ . Inserting (30) into (31) to (33), we can then attain the following PMI-LCP contact formulation:

$$\begin{array}{c}\left(\begin{array}{c}{\sigma }_k^n\hfill \\ {\sigma }_k^t\hfill \\ {\sigma }_k^f\hfill \end{array}\right)=\left[\begin{array}{ccc}{\bar{N}}_k^T{A}_k^{-1}{\bar{N}}_k\hfill & {\bar{N}}_k^T{A}_k^{-1}{\bar{D}}_k\hfill & 0\hfill \\ {\bar{D}}_k^T{A}_k^{-1}{\bar{N}}_k\hfill & {\bar{D}}_k^T{A}_k^{-1}{\bar{D}}_k\hfill & U\hfill \\ \mu I_{l_c\times l_c}\hfill & -U^T\hfill & 0\hfill \end{array}\right]\left(\begin{array}{c}{\lambda }_k^n\hfill \\ {\lambda }_k^t\hfill \\ {\lambda }_k^f\hfill \end{array}\right)+\left(\begin{array}{c}{\gamma }_k^n\hfill \\ {\gamma }_k^t\hfill \\ 0\hfill \end{array}\right)\\ {\sigma }_k^n,{\sigma }_k^t,{\sigma }_k^f,{\lambda }_k^n,{\lambda }_k^t,{\lambda }_k^f\ge 0\\ {\sigma }_k^{\mathrm{nT}}{\lambda }_k^n={\sigma }_k^{\mathrm{tT}}{\lambda }_k^t={\sigma }_k^{\mathrm{fT}}{\lambda }_k^f=0\end{array}$$

(34)

where ${\gamma }_k^n:={\bar{N}}_k^T{A}_k^{-1}{B}_k$ and ${\gamma }_k^t={\bar{D}}_k^T{A}_k^{-1}{B}_k$ . This contact PMI-LCP (34) can be solved by using the interior point method or the Lemke pivoting algorithm (Cottle et al., 2009).

The complementarity conditions (31) to (33) naturally suggest passivity of our contact PMI-LCP formulation (34), that is, it is lossless if no tangential slip occurs, or dissipative if it does occur. This passivity of PMI-LCP (34) is clearly exhibited in Figure 6, where a rigid box free-falls on a rigid level surface and bounces back and forth, with the total energy being constant during the flight stages, but strictly decreasing during the contact stages, during which dissipative tangential slips occur. Note that this simulates lossless (i.e. elastic) collision along the normal direction and dissipative slips along the tangential directions, with no deformation allowed between the two rigid objects. With no slips, this would result in completely lossless/elastic collision.

Fig. 6.

Profiles of total energy and height with PMI-LCP (34): elastic–normal and dissipative–tangential slip contact attained with total energy remaining constant during flight, but strictly decreased during contact. Sudden changes of H around 4 s to 4.2 s are due to collision with flipping.

This ability of our PMI-LCP (34) to naturally render (marginally stable) elastic contact is in stark contrast to other contact rendering results based on IEI or SEI, which typically can only produce plastic contact unless some extra velocity restitution models are deployed. This is in fact due to the adoption of $v_{k+1}$ as the representative velocity by the IEI and SEI (see Section 2.3). See Figure 7, where the same box-falling scenario as in Figure 6 results in purely plastic collision with the SEI. Of course, in many applications plastic collisions are necessary, which can also be easily achieved by our PMI-LCP (34) simply by replacing ${\widehat{V}}_k^a$ in ${\sigma }_k^n,{\sigma }_k^t$ by $V_{k+1}$ (see Figure 8).

Fig. 7.

Profiles of total energy and height with SEI: although the setting is the same as in Figure 6, completely plastic collision is attained (see also Figure 5).

Fig. 8.

Profiles of total energy and height under PMI-LCP (34) with ${\widehat{V}}_k^a$ replaced by $V_{k+1}$ for (31) and (32): completely plastic collision attained as in Figure 7 under SEI.

Here, note that, in contrast to Section 3.2 (i.e. approximate passivity of ${\lambda }_k$ (25)), this passivity of the PMI-LCP is exact, since we do not adopt the Jacobian approximation used in (22). Consequently, the damping injection, added to enforce passivity of the approximate potential action in Section 3.2, is not necessary here either. Our PMI-LCP algorithm is formulated at the velocity level (i.e. (27)), and thus suffers from the issue of constraint drift (induced by initial collision detection error or numerical precision error). This constraint drift can be easily remedied by incorporating a drift-correction term into ${\gamma }_k^n$ in (34) as done in a similar manner in Otaduy et al. (2009). Another interesting direction for this is to combine the PMI-LCP formulation with the approximate action of the potential ${\psi }_{{\lambda }_k}$ of Section 3.2 to stabilize this constraint drift. Passivity (and, thus, stability) of these drift corrections or combination of the PMI-LCP with its approximate constraint potential, however, has not been established and is a topic for future research.

4. PMI-based simulation in generalized coordinates

In this section, we derive the PMI of the Lagrange dynamics (4) in generalized coordinates. We first show that the PMI formulations developed in Section 2.2 and Section 3.1 are not able to enforce passivity for the Lagrange dynamics (4) due to the nonlinearity of the dynamics. To overcome this, we propose a certain coordinate transformation, with which the nonlinear Lagrange dynamics (4) can be converted into a form compatible with the PMI formalism to enforce passivity in the discrete-time domain. We further show how to derive the PMI-LCP formulation in generalized coordinates for the Lagrange dynamics (4) to simulate point Coulomb friction contact, as done in Section 3.3 for the PMI in maximal coordinates.

4.1. PMI of Lagrange dynamics

Consider the multi-DOF nonlinear Lagrange dynamics (4) with the generalized coordinates $q\in {\Re }^n$ . Here, we focus on quadratic potential functions $\psi (q)=\frac{1}{2}{q}^T\mathrm{Kq}$ , since with them we can establish discrete-time passivity exactly with no approximations necessary. Of course, other types of analytic nonlinear potentials $\psi (q)$ can also be used in the same manner, with their potential action approximated as in Section 3.2 to preserve the non-iterativeness of the PMI simulation and with some damping B injection to enforce discrete-time passivity even with this approximation. Contact with external objects, whose positions can be nonlinear with respect to q (i.e. on-link contact), is to be rendered by the PMI-LCP formulation of Section 4.3 via $f_c$ in (4) instead of via the potential $\psi (q)$ .

Now, suppose M(q) in (4) be constant (i.e. $M(q)=M$ ). This means $C(q,\stackrel{.}{q})$ in (4) is skew-symmetric, since $Ṁ-2C$ is skew-symmetric. We can then write the PMI expression of (4) as a vector version of (9), with $x_k,v_k$ replaced by $q_k,v_k:={\stackrel{.}{q}}_k\in {\Re }^n$ and also with $C_k{\widehat{v}}_k:=C(q_k,v_k){\widehat{v}}_k$ added in (9) similarly to the term $S(J{w}_k){ŵ}_k$ in (13). Further, similar to (15), we can show that this PMI of Lagrange dynamics (4) with constant M and linear spring potential $d{\psi }_k^T=K{q}_k$ is discrete-time passive, that is, $\forall N\ge 1$ ,

$$\begin{array}{c}\sum _{k=0}^{N-1}{[f_k+f_k^c]}^T{\widehat{v}}_k{T}_k=\sum _{k=0}^{N-1}\left[E_{k+1}-E_k+{\widehat{v}}_k^T{C}_k{\widehat{v}}_k{T}_k\right]\hfill \\ =E_{\bar{N}+1}-E_o\ge -E_o\end{array}$$

(35)

where $C_k$ is skew-symmetric, ${\widehat{v}}_k:=(v_{k+1}+v_k)∕2$ is the “representative” velocity during $T_k$ , and $E_k:=\frac{1}{2}{v}_k^{T}M{v}_k+\frac{1}{2}{q}_k^{T}K{q}_k$ is the total energy.

This discrete-time passivity (35), however, is not granted for non-constant M. That is, from M(q) being configuration-dependent, we have $M_{k+1}:=M(q_{k+1})\ne M_k:=M(q_k)$ , with which, instead of (35), we have

$$\begin{array}{c}\sum _{k=0}^{N-1}[f_k^c+f_k{]}^T{\widehat{v}}_k{T}_k=\frac{1}{2}{q}_N^{T}K{q}_N-\frac{1}{2}{q}_o^{T}K{q}_o\hfill \\ +\sum _{k=0}^{N-1}\left[\frac{v_{k+1}^T{M}_k{v}_{k+1}}{2}-\frac{v_k^T{M}_k{v}_k}{2}+{\widehat{v}}_k^T{C}_k{\widehat{v}}_k{T}_k\right]\end{array}$$

where, differently from (35), the kinetic energies in the brackets do not cancel each other out over the successive $T_k$ and $T_{k+1}$ , and $C_k$ is not skew-symmetric either. This is because the skew-symmetry of $Ṁ-2C$ (or, equivalently, $Ṁ=C+C^T$ ) cannot be simulated during $T_k$ , since doing so requires future information $M_{k+1},C_{k+1}$ at the onset of $T_k$ . To address this, here we utilize the coordinate transformation proposed in Lee et al. (2012), with which the nonlinear Lagrange dynamics (4) is transformed intoanother Lagrangian form with constant inertia and a skew-symmetric Coriolis term, and thus can still be passive by using PMI similarly to (13).

More precisely, define the following (continuous-time) coordinate transformation:

$$\begin{array}{c}\xi :=M^{\frac{1}{2}}(q)\stackrel{.}{q},\delta :=M^{-\frac{1}{2}}(q)f\hfill \end{array}$$

(36)

with which the kinetic energy can be written as

$$\frac{1}{2}{\stackrel{.}{q}}^{T}M(q)\stackrel{.}{q}=\frac{1}{2}{\xi }^T\mathrm{I\xi }$$

where $I\in {\Re }^{n\times n}$ is the identity matrix, implying that, for the transformed velocity $\xi$ , the inertia will be this constant I. Using

$$I\stackrel{.}{\xi }=M^{\frac{1}{2}}(q)\stackrel{··}{q}+\frac{d}{\mathrm{dt}}[M^{\frac{1}{2}}(q)]\stackrel{.}{q}$$

we can then rewrite (4):

$$\begin{array}{c}I\stackrel{.}{\xi }+Q(q,\stackrel{.}{q})\xi +M^{-\frac{1}{2}}(q)\mathrm{Kq}={\delta }_c+\delta \hfill \end{array}$$

(37)

where the Coriolis-like term $Q(q,\stackrel{.}{q})\in {\Re }^{n\times n}$ is defined by

$$\begin{array}{c}Q(q,\stackrel{.}{q}):=M^{-\frac{1}{2}}C{M}^{-\frac{1}{2}}-\frac{d}{\mathrm{dt}}[M^{\frac{1}{2}}]M^{-\frac{1}{2}}\hfill \end{array}$$

(38)

with arguments omitted for brevity.

This matrix Q in (38) is in fact skew-symmetric. That is, with M being symmetric, we have

$$\begin{array}{c}Q+Q^T=M^{-\frac{1}{2}}[C+C^T]M^{-\frac{1}{2}}-\frac{d}{\mathrm{dt}}[M^{\frac{1}{2}}]M^{-\frac{1}{2}}-M^{-\frac{1}{2}}\frac{d}{\mathrm{dt}}[M^{\frac{1}{2}}]\hfill \\ =M^{-\frac{1}{2}}\frac{d}{\mathrm{dt}}[M^{\frac{1}{2}}{M}^{\frac{1}{2}}]M^{-\frac{1}{2}}-\frac{d}{\mathrm{dt}}[M^{\frac{1}{2}}]M^{-\frac{1}{2}}-M^{-\frac{1}{2}}\frac{d}{\mathrm{dt}}[M^{\frac{1}{2}}]\\ =0\end{array}$$

where we use the fact that $Ṁ=C+C^T$ from $Ṁ-2C$ being skew-symmetric. We can also show that the transformed dynamics (37) preserves the same passivity as (5): $\forall T\ge 0$ ,

$$\begin{array}{c}{\int }_0^T{[f_c+f]}^T\stackrel{.}{q}\mathrm{dt}={\int }_0^T{[{\delta }_c+\delta ]}^T\xi \mathrm{dt}\hfill \\ ={\int }_0^T\left[{\xi }^{T}I\stackrel{.}{\xi }+{\xi }^T\mathrm{Q\xi }+{\xi }^T{M}^{-\frac{1}{2}}\mathrm{Kq}\right]\mathrm{dt}\\ =E(T)-E(0)\ge -E(0)\end{array}$$

where $E:=\frac{1}{2}{\stackrel{.}{q}}^{T}M(q)\stackrel{.}{q}+\frac{1}{2}{q}^T\mathrm{Kq}=\frac{1}{2}{\xi }^T\mathrm{I\xi }+\frac{1}{2}{q}^T\mathrm{Kq}$ is the total energy.

Since the transformed dynamics (37) has constant inertia I and skew-symmetric Q, we can now passively simulate (37) similarly to (9) and (13):

$$\begin{array}{c}I\frac{{\xi }_{k+1}-{\xi }_k}{T_k}+Q_k\frac{{\xi }_{k+1}+{\xi }_k}{2}+M_k^{-\frac{1}{2}}K\frac{q_{k+1}+q_k}{2}\hfill \\ ={\delta }_{c_k}+{\delta }_k\end{array}$$

(39)

where ${\delta }_{c_k}:=M_k^{-\frac{1}{2}}{f}_{c_k}$ , ${\delta }_k:=M_k^{-\frac{1}{2}}{f}_k$ , and $Q_k:=Q(q_k,v_k)$ . We can then show that

$$\begin{array}{c}\sum _{k=1}^{N-1}{[{\delta }_{c_k}+{\delta }_k]}^T{\widehat{\xi }}_k{T}_k\\ =\sum _{k=0}^{N-1}\left[{\kappa }_{k+1}-{\kappa }_k+{\widehat{\xi }}_k^T{Q}_k{\widehat{\xi }}_k{T}_k+{\widehat{\xi }}_k^T{M}_k^{-\frac{1}{2}}K\frac{q_{k+1}+q_k}{2}{T}_k\right]\\ ={\kappa }_N-{\kappa }_o+\sum _{k=0}^{N-1}{\widehat{\xi }}_k^T{M}_k^{-\frac{1}{2}}K\frac{q_{k+1}+q_k}{2}{T}_k\end{array}$$

where ${\widehat{\xi }}_k:=({\xi }_{k+1}+{\xi }_k)∕2$ , ${\kappa }_k:={\xi }_k^{T}I{\xi }_k∕2=v_k^T{M}_k{v}_k∕2$ , and we use the fact that $Q_k$ is skew-symmetric. This then also suggests the following kinematics law, similar to (8), but transformed via $M_k^{\frac{1}{2}}$ :

$$\begin{array}{c}\frac{q_{k+1}-q_k}{T_k}=M_k^{-\frac{1}{2}}\frac{{\xi }_{k+1}+{\xi }_k}{2}\hfill \end{array}$$

(40)

which is indeed consistent with the continuous-time transformation (36).

Note that, given the continuous-time Lagrangian dynamics (4), we first transform it to (37), then discretize it via (39), and then connect to the configuration q of the original dynamics through (40). To compute $Q_k$ via (38), we may also use the following (discrete-time) Sylvester equation (Horn and Johnson, 2012) to obtain ${Ṁ}_k^{\frac{1}{2}}$ . That is, given $M_k^{\frac{1}{2}},C_k$ , find ${Ṁ}_k^{\frac{1}{2}}$ such that

$${Ṁ}_k^{\frac{1}{2}}{M}_k^{\frac{1}{2}}+M_k^{\frac{1}{2}}{Ṁ}_k^{\frac{1}{2}}={Ṁ}_k=[C_k+C_k^T]$$

with which we can also show that $Q_k=M_k^{-\frac{1}{2}}{C}_k{M}_k^{-\frac{1}{2}}-{Ṁ}_k^{\frac{1}{2}}{M}_k^{-\frac{1}{2}}$ is skew-symmetric. In Section 4.3, we will explain how to render the contact force ${\delta }_{c_k}$ in (39) based on the point contact Coulomb friction model by incorporating the LCP formulation into the PMI expression (39) as similarly done in Section 3.3. Before doing so, let us first compare the PMI with symplectic integrators to shed light on some unique and useful properties of the PMI.

4.2. Comparison of PMI with symplectic integrators

To define symplectic integrators for the Lagrangian dynamics (4), let us define its Hamiltonian,

$$\begin{array}{c}H(p,q)=p^{T}M{(q)}^{-1}p+U(q)\hfill \end{array}$$

where $q\in {\Re }^n$ is the configuration, $p:=M(q)\stackrel{.}{q}\in {\Re }^n$ is the conjugate momentum, and $U(q)\in \Re$ is the potential energy. Then, the Hamiltonian dynamics formulation of (4) is given by

$$\begin{array}{c}ṗ=-\frac{\partial H}{\partial q}(p,q),\stackrel{.}{q}=\frac{\partial H}{\partial p}(p,q)\hfill \end{array}$$

which can be written in the following matrix form:

$$\begin{array}{c}ẏ:=\left(\begin{array}{c}ṗ\hfill \\ \stackrel{°}{q}\hfill \end{array}\right)={\left[\begin{array}{cc}0\hfill & I\hfill \\ -I\hfill & 0\hfill \end{array}\right]}^{-1}\left(\begin{array}{c}\frac{\partial H}{\partial p}\hfill \\ \frac{\partial H}{\partial q}\hfill \end{array}\right)=:{\mathsf{J}}^{-1}\cdot \nabla H(y)\hfill \end{array}$$

(41)

where $y:=[p;q]\in {\Re }^{2n}$ is the coordinates in the tangent bundle, $\nabla H(y)=[\frac{\partial H}{\partial p};\frac{\partial H}{\partial q}]\in {\Re }^{2n}$ is the gradient of H, and $\mathsf{J}\in {\Re }^{2n\times 2n}$ is the (skew-symmetric) canonical structure matrix. The solution of the Hamiltonian system (41) can then be written using its flow map ${\psi }_t$ with ${\psi }_t(y_0):=y(t,y_0)$ .

The Hamiltonian system (41) possesses three fundamental properties, namely, energy conservation (or passivity), momentum conservation, and symplecticity. Since the first two are familiar, here we explain only the symplecticity. A linear mapping $A:{\Re }^{2n}\to {\Re }^{2n}$ is called symplectic if

$$\begin{array}{c}w(A{\zeta }_1,A{\zeta }_2)=w({\zeta }_1,{\zeta }_2)\hfill \end{array}$$

(42)

for ${\zeta }_1:=[{\zeta }_1^p,{\zeta }_1^q],{\zeta }_2:=[{\zeta }_2^p,{\zeta }_2^q]\in {\Re }^{2n}$ with ${\zeta }_1^p,{\zeta }_1^q,{\zeta }_2^p,{\zeta }_2^q\in {\Re }^n$ , where

$$\begin{array}{c}w({\zeta }_1,{\zeta }_2):=\sum _{i=1}^{n}det\left(\begin{array}{cc}{\zeta }_{1,i}^p\hfill & {\zeta }_{2,i}^p\hfill \\ {\zeta }_{1,i}^q\hfill & {\zeta }_{2,i}^q\hfill \end{array}\right)=\sum _{i=1}^n({\zeta }_{1,i}^p{\zeta }_{2,i}^q-{\zeta }_{1,i}^q{\zeta }_{2,i}^p)\hfill \end{array}$$

which is the sum of the oriented areas. In other words, symplecticity implies area preservation through the map A. For the Hamiltonian system (41), the symplecticity condition is given by

$$\begin{array}{c}{\left(\frac{\partial {\psi }_t}{\partial y_0}\right)}^T\mathsf{J}\left(\frac{\partial {\psi }_t}{\partial y_0}\right)=\mathsf{J}\hfill \end{array}$$

Symplectic integrators are those that preserve this symplecticity property of the continuous-time Hamiltonian system (41). The following symplectic Euler integration and implicit midpoint integration are those symplectic integrators:

$$\begin{array}{cc}{y}_{k+1}\hfill & =y_k+{\mathsf{J}}^{-1}\cdot \nabla H(p_{k+1},q_k)T\\ y_{k+1}\hfill & =y_k+{\mathsf{J}}^{-1}\cdot \nabla H(p_k,q_{k+1})T\end{array}$$

(43)

$$\begin{array}{cc}{y}_{k+1}\hfill & =y_k+{\mathsf{J}}^{-1}\cdot \nabla H({\widehat{p}}_k,{\widehat{q}}_k)T\end{array}$$

(44)

with the first two being SEI and the last being IMI, where $y_k:=[p_k;q_k]\in {\Re }^{2n}$ , $T>0$ is the integration step, and ${\widehat{\star }}_k:=\frac{{\star }_{k+1}+{\star }_k}{2}$ . Both of these SEI and IMI are known to be symplectic and momentum-conserving, but not energy-preserving, where the symplecticity condition is given by

$$\begin{array}{c}{\left(\frac{\partial y_{k+1}}{\partial y_k}\right)}^T\mathsf{J}\left(\frac{\partial y_{k+1}}{\partial y_k}\right)=\mathsf{J}\hfill \end{array}$$

(45)

Now, let us see if the PMI (39)–(40) is also symplectic. For this, from (39) and (40), we can obtain the following equation:

$$\begin{array}{c}\left(\begin{array}{c}{p}_{k+1}\hfill \\ q_{k+1}\hfill \end{array}\right)=\left[\begin{array}{cc}{M}_{k+1}^{\frac{1}{2}}(q_{k+1})\hfill & 0\hfill \\ 0\hfill & I\hfill \end{array}\right]{\left(\frac{I}{T_k}+S_k(p_k,q_k)\right)}^{-1}\hfill \\ \times \left(\frac{I}{T_k}-S_k(p_k,q_k)\right)\left[\begin{array}{cc}{M}_k^{-\frac{1}{2}}(q_k)\hfill & 0\hfill \\ 0\hfill & I\hfill \end{array}\right]\left(\begin{array}{c}{p}_k\hfill \\ q_k\hfill \end{array}\right)\end{array}$$

(46)

where $p_k=M^{\frac{1}{2}}(q_k){\xi }_k\in {\Re }^n$ and $S_k(p_k,q_k):=\frac{1}{2}\left[Q_k(p_k,q_k),M_k^{-\frac{1}{2}}(q_k)K;M_k^{-\frac{1}{2}}(q_k),0\right]\in {\Re }^{2n\times 2n}$ .

To derive (45) from (46) turns out to be too complicated to analytically perform. Note however that we have $M_{k+1}^{\frac{1}{2}}$ in (46), which in general can take any positive-definite matrix, suggesting the PMI would be unlikely to be symplectic. For this, here we resort to a numerical approach, that is, we measure how the area of a certain shape in the space of $(p,q)\in {\Re }^{2n}$ evolves in the course of the PMI integration of (39) and (40). The results are shown in Figure 9, where we can see that the PMI is energy-preserving (from passivity), but not symplectic, although it still stays near being symplectic. We can show that the PMI is not momentum-preserving either. This then means that the PMI is not the IMI, although for linear mechanical systems they become the same.

Fig. 9.

Simulation of two-serial-link robot arm with linear joint springs with SEI and PMI: the SEI is symplectic, but not exactly energy-preserving, whereas the PMI is energy-preserving, but not exactly symplectic.

Even so, we believe the PMI has advantages for haptic rendering and interactive simulation as compared to the symplectic integrators, particularly due to the following two powerful properties: (1) non-iterativeness and (2) passivity. Note first from (43) and (44) that the symplectic integrators should in general solve the nonlinear map (i.e. $\nabla H$ ) in real time in an iterative manner due to the presence of the same variables both on the left-hand side and inside of $\nabla H(\star )$ in (43) and (44). This iterative nature of the symplectic integrators, combined with their non-passivity (i.e. instability possible), poses a challenging issue in practice, that is, this nonlinear map must be solved precisely to maintain stability, which may be time-consuming for complex $\nabla H$ or require only a small time-step T to ensure sufficiently fastconvergence of this nonlinear map solving, which are both not very suitable for real-time simulation interaction. From its being dissipative, the IEI, in contrast, is known to be less sensitive to this solution precision and is often chosen over SEI or IMI when this issue becomes important.

On the other hand, note from (26) or (39) and (40) that the PMI is non-iterative with no requirement to solve such a nonlinear map. This non-iterativeness of the PMI is attained due to the coordinate transformation (36) (for nonlinearity from inertia) and the potential action approximation (e.g. (24) for potential nonlinearity). Of course, with this approximate action of nonlinear potentials, the passivity will also in general be only approximate. This “leakage” of passivity can however be rather easily resolved by injecting (typically small) linear damping B, which can still be simulated quickly due to the non-iterativeness of the PMI and whose passivity is also still guaranteed due to the passivity of the PMI. This unconditionally stable/passive damping allows for the so-called “controllable damping” (e.g. increase (or respectively decrease) B for more stable (or respectively less stable) behavior while preserving simulation stability), which is not possible for the symplectic integrators, theoretically or practically (i.e. precision of nonlinear map-solving). The non-iterativeness and passivity, we believe, are more relevant than symplecticity and momentum-preserving to haptic rendering and interactive simulation, since stability and real-time interactivity are their most fundamental requirements and their simulated objects are also typically under the influence of frequent human or contact forces, which is different from objects in celestial mechanics or molecular simulation, where long-term isolated behavior simulation is more important.

4.3. PMI-LCP for multi-point Coulomb friction contact

Following the procedure in Section 3.3, here we utilize the point contact model with Coulomb friction. Here, we assume the multi-DOF Lagrange system (39) interacts with a static virtual rigid object. We choose this static object for simplicity: the derivation presented below can be easily applied to multiple moving rigid objects, which can then be rendered by using the PMI expression in maximal coordinates from Section 3. We also adopt the velocity-level complementary condition as in (27) and the $\bar{p}$ -pyramidal cone approximation as in (29), and denote by $l_c$ the number of contact points between the Lagrange system (39) and the object.

Then, the contact force ${\delta }_{i,k}^c\in {\Re }^n$ at the i th contact point ( $i\in \{1,2,\dots ,l_c\}$ ) can be written as

$$\begin{array}{c}{\delta }_{i,k}^c=M_k^{-\frac{1}{2}}{f}_{i,k}^c=M_k^{-\frac{1}{2}}{J}_{c_i,k}^T\left[N_{i,k}{\lambda }_{i,k}^n+D_{i,k}{\lambda }_{i,k}^t\right]\hfill \end{array}$$

(47)

where $J_{c_i,k}(q_k)\in {\Re }^{3\times n}$ is the Jacobian of the i th contact point, and ${\lambda }_{i,k}^n\in \Re$ and ${\lambda }_{i,k}^t\in {\Re }^{\bar{p}}$ are the normal contact and tangential friction forces with $N_{i,k}\in {\Re }^3$ and $D_{i,k}\in {\Re }^{3\times \bar{p}}$ respectively encoding their directions in ${\Re }^3$ , defined in a manner similar to (29). Combining (47) for all $l_c$ contact points, we can then obtain the total contact force ${\delta }_{c_k}$ in (39):

$${\delta }_{c_k}={\bar{N}}_k{\lambda }_k^n+{\bar{D}}_k{\lambda }_k^t$$

where

$$\begin{array}{c}{\lambda }_k^n:=[{\lambda }_{1,k}^n;{\lambda }_{2,k}^n;\dots ;{\lambda }_{l_c,k}^n]\in {\Re }^{l_c}\hfill \\ {\lambda }_k^t:=[{\lambda }_{1,k}^t;{\lambda }_{2,k}^t;\dots ;{\lambda }_{l_c,k}^t]\in {\Re }^{l_c\times \bar{p}}\hfill \\ {\bar{N}}_k:=M_k^{-\frac{1}{2}}\left[J_{c_1,k}^T{N}_{1,k},J_{c_2,k}^T{N}_{2,k},\dots J_{c_{l_c},k}^T{N}_{l_c,k}\right]\in {\Re }^{n\times l_c}\hfill \\ {\bar{D}}_k:=M_k^{-\frac{1}{2}}\left[J_{c_1,k}^T{D}_{1,k},J_{c_2,k}^T{D}_{2,k},\dots J_{c_{l_c},k}^T{D}_{l_c,k}\right]\in {\Re }^{n\times (l_c\times \bar{p})}\hfill \end{array}$$

with which we can also rewrite the PMI expression of the Lagrange dynamics (39) similarly to (30), that is,

$$\begin{array}{cc}{\widehat{\xi }}_k\hfill & =A_k^{-1}\left[B_k+{\bar{N}}_k{\lambda }_k^n+{\bar{D}}_k{\lambda }_k^t\right]\end{array}$$

(48)

where $A_k:=\frac{2}{T_k}I+Q_k+\frac{1}{2}{M}_k^{-\frac{1}{2}}K{M}_k^{-\frac{1}{2}}{T}_k$ and $B_k:=\frac{2}{T_k}{\xi }_k-M_k^{-\frac{1}{2}}K{q}_k$ .

Given this ${\widehat{\xi }}_k$ in (48), we can compute the velocity of all $l_c$ contact points:

$$\begin{array}{c}{\widehat{V}}_k^c=J_{c_k}{M}_k^{-\frac{1}{2}}{\widehat{\xi }}_k\in {\Re }^{3l_c}\hfill \end{array}$$

(49)

where $J_{c_k}:=\left[J_{c_1,k};J_{c_2,k};\dots J_{c_{l_c},k}\right]\in {\Re }^{3l_c\times n}$ is the combined Jacobian of the $l_c$ contact points. The contact complementarity conditions (Lloyd, 2005) can be written similarly to (31) to (33):

$$\begin{array}{c}0\le {\sigma }_k^n\perp {\lambda }_k^n\ge 0\\ 0\le {\sigma }_k^t\perp {\lambda }_k^t\ge 0\\ 0\le {\sigma }_k^f\perp {\lambda }_k^f\ge 0\end{array}$$

where

$$\begin{array}{c}{\sigma }_k^n:=N_k^T{\widehat{V}}_k^c\in {\Re }^{l_c}\\ {\sigma }_k^t:=D_k^T{\widehat{V}}_k^c+U{\lambda }_k^f\in {\Re }^{l_c\times \bar{p}}\\ {\sigma }_k^f:=\bar{\mu }{\lambda }_k^n-U^T{\lambda }_k^t\in {\Re }^{l_c}\end{array}$$

where $N_k:=\mathrm{tr}[N_{1,k},N_{2,k},\dots N_{l_c,k}]\in {\Re }^{3l_c\times l_c}$ , $D_k:=\mathrm{tr}[D_{1,k},D_{2,k},\dots D_{l_c,k}]\in {\Re }^{3l_c\times (l_c\times \bar{p})}$ , and $\bar{\mu }:=\mu I_{l_c\times l_c}$ and $U:=\mathrm{tr}[1_{\bar{p}},1_{\bar{p}},..1_{\bar{p}}]\in {\Re }^{(l_c\times \bar{p})\times \bar{p}}$ as for (31) to (33).

Finally, similarly to (34), we can construct the PMI-LCP formulation in generalized coordinates:

$$\begin{array}{c}\left(\begin{array}{c}{\sigma }_k^n\hfill \\ {\sigma }_k^t\hfill \\ {\sigma }_k^f\hfill \end{array}\right)=\left[\begin{array}{ccc}{W}_{11}\hfill & W_{12}\hfill & 0\hfill \\ W_{21}\hfill & W_{22}\hfill & U\hfill \\ \mu I_{l_c\times l_c}\hfill & -U^T\hfill & 0\hfill \end{array}\right]\left(\begin{array}{c}{\lambda }_k^n\hfill \\ {\lambda }_k^t\hfill \\ {\lambda }_k^f\hfill \end{array}\right)+\left(\begin{array}{c}{\gamma }_k^n\hfill \\ {\gamma }_k^t\hfill \\ 0\hfill \end{array}\right)\\ {\sigma }_k^n,{\sigma }_k^t,{\sigma }_k^f,{\lambda }_k^n,{\lambda }_k^t,{\lambda }_k^f\ge 0\\ {\sigma }_k^{\mathrm{nT}}{\lambda }_k^n={\sigma }_k^{\mathrm{tT}}{\lambda }_k^t={\sigma }_k^{\mathrm{fT}}{\lambda }_k^f=0\end{array}$$

(50)

where ${\gamma }_k^n:=N_k^T{J}_{c_k}{M}_k^{-\frac{1}{2}}{A}_k^{-1}{B}_k$ , ${\gamma }_k^t:=D_k^T{J}_{c_k}{M}_k^{-\frac{1}{2}}{A}_k^{-1}{B}_k$ and

$$\begin{array}{c}{W}_{11}:=N_k^T{J}_{c_k}{M}_k^{-\frac{1}{2}}{A}_k^{-1}{M}_k^{-\frac{1}{2}}{J}_{c_k}^T{N}_k\hfill \\ W_{12}:=N_k^T{J}_{c_k}{M}_k^{-\frac{1}{2}}{A}_k^{-1}{M}_k^{-\frac{1}{2}}{J}_{c_k}^T{D}_k=W_{21}^T\hfill \\ W_{22}:=D_k^T{J}_{c_k}{M}_k^{-\frac{1}{2}}{A}_k^{-1}{M}_k^{-\frac{1}{2}}{J}_{c_k}^T{D}_k\hfill \end{array}$$

which can be solved by an interior point method or Lemke’s pivoting algorithm (Cottle et al., 2009).

5. Illustrative examples

So far, we have derived the PMI-based simulation framework, which, due to its enforcing passivity being effective in practice, can stably simulate a wide range of parameters (e.g. very light/stiff objects with no damping) regardless of (possibly varying/slow) integration steps $T_k$ , which can occur during, for example, intermittent multi-point collisions among virtual objects. To validate and demonstrate these (and other) advantages of our proposed PMI-based simulation framework, here we perform three illustrative examples, namely haptic operation of a peg-in-hole task (Section 5.1), haptic interaction with a flexible beam (Section 5.2), and multi-point power grasping of a multi-fingered UATD hand (Section 5.3).

5.1. Peg-in-hole task

For the peg-in-hole task, we model the peg and the box with holes as rigid bodies, and the two multi-fingered virtual hands as articulated rigid bodies with joint constraints as shown in Figure 4, each composed of 10 articulated rigid bodies. To simulate these, we apply the PMI-based formulation in maximal coordinates from Section 3. We also utilize the wearable cutaneous haptic device (CHD) proposed in Kim et al. (2016), which can provide 3-DOF contact force feedback at the fingertip. This CHD is different from typical haptic devices (e.g. Omega 3, Geomagic Phantom), as it cannot generate kinesthetic force feedback. Consequently, we utilize unilateral virtual coupling to allow the user to command the position/orientation of the virtual hands/fingers, whereas the fingertip contact force is directly fed back to the CHD.

More precisely, the SE(3) springs, consisting of the E(3) linear spring in Section 2.2 and the SO(3) spring in Section 3.1, are used as uni-directional position/orientation-commanding virtual coupling (e.g. Kim and Lee, 2017; Lee, 2009). On the other hand, the contact force of each fingertip of the virtual hands during the manipulation (computed via the PMI-LCP in Section 3.3) is directly fed back to the cutaneous fingertip haptic device. For this, we use

$$\begin{array}{c}{f}_{f,k}=R_{f,k}^T{H}_{\mathrm{filter}}(f_{i,k-1}^c,f_{i,k-2}^c,\dots f_{i,k-p_f}^c)\hfill \end{array}$$

where $R_{f,k}$ is the rotation matrix of the fingertip body frame with respect to the inertial frame, and $H_{\mathrm{filter}}$ is a certain low-pass filter to process the $p_f$ -number of the computed contact force data $f_{i,k-j}^c$ , $j=1,\dots ,p_f$ . The filter $H_{\mathrm{filter}}$ is used here to provide smooth cutaneous haptic feedback $f_{f,k}$ (with 1 kHz) from the computed contact force $f_{i,k}^c$ , whose update rate varies between 30 Hz and 300 Hz due to the simulation of the large number of rigid bodies (i.e. two hands, each composed of 10 rigid bodies, one box with holes, and one peg) and multi-point contacts between them. Note that we cannot use $f_{i,k}^c$ for $H_{\mathrm{filter}}$ , since it can be determined only after $T_k=[t_k,t_{k+1})$ . The computational lag from $H_{\mathrm{filter}}$ and the unilateralness of the virtual coupling do not affect stability, since the device–human coupling here is only cutaneous, not kinesthetic (Prattichizzo et al., 2013).

Snapshots of this peg-in-hole task are shown in Figure 10 with its fingertip forces also presented in Figure 11. We set the mass and inertia of the rigid bodies of the hand as $m=0.001$ kg and $J=m{I}_{3\times 3}$ kgm², and the stiffness of the spring as $k_{\lambda }=100$ N/m, with the stiffness/mass ratio having a fairly high value of $10^5$ . We found that, with these parameters, the cutaneous haptic feedback rendering is transparent enough, since (1) the virtual hand can (unilaterally) follow the user command fast enough with this high stiffness/mass ratio of $10^5$ , and (2) the cutaneous force feedback $f_{f,k}$ , as computed above, is directly displayed to the fingertip with no intermediate compliance. Note that this differs from the case of kinesthetic devices with bilateral virtual coupling, where, for the transparency, the virtual coupling should typically be much higher, as it defines bilateral compliant coupling between the virtual proxy and the real haptic device. Of course, we can also adopt kinesthetic haptic devices here along with, for example, the virtual coupling developed for the PMI (Huang and Lee, 2011; Lee et al., 2012), or passive set-position modulation (Lee and Huang, 2010). In this case, even if we use the PMI to enforce the passivity of the virtual worlds, the (hybrid/bilateral) virtual coupling gain will be still limited by the device’s physical damping; see Huang and Lee (2011).

Fig. 10.

Snapshots of peg-in-hole task with two multi-fingered hands, each modeled as an articulated rigid body as shown in Figure 4 and making contact with the objects via friction contact simulated by the PMI-LCP in Section 3.3.

Fig. 11.

Fingertip contact force during the peg-in-hole task of Figure 10 expressed in the fingertip coordinate frame (with the z-axis being the contact normal direction).

5.2. Haptic interaction with flexible beam

We model the flexible beam to be an aluminum beam with the dimensions $(0.907,0.05,0.012)$ m and a mass of 1.4683 kg. We also model its dynamics by using nodal coordinates with the three natural frequencies $(11.91,74.65,209.03)$ Hz. Here, for simplicity, we assume the beam motion to be planar. Then, the beam dynamics can be modeled in mixed maximal–generalized coordinates and by using an Euler-Bernoulli model (Rao and Yap, 1995) such that

$$\begin{array}{c}\left[\begin{array}{cc}{M}_{\mathrm{xx}}\hfill & M_{x\delta }\hfill \\ M_{\delta x}\hfill & M_{\mathrm{\delta \delta }}\hfill \end{array}\right]\left(\begin{array}{c}ẍ\hfill \\ \stackrel{··}{\delta }\hfill \end{array}\right)+\left[\begin{array}{cc}{C}_{\mathrm{xx}}\hfill & C_{x\delta }\hfill \\ C_{\delta x}\hfill & 0\hfill \end{array}\right]\left(\begin{array}{c}ẋ\hfill \\ \stackrel{°}{\delta }\hfill \end{array}\right)+d{\psi }^T(\delta )=f\hfill \end{array}$$

where $x\in {\Re }^3$ is the maximal coordinate of the beam end in SE(2), $\delta \in {\Re }^3$ is the generalized coordinate representing each vibration mode, $d{\psi }^T(\delta ):=K_{q}q:=[0;K_{\mathrm{\delta \delta }}\delta ]$ , and $M_{\mathrm{\delta \delta }},K_{\mathrm{\delta \delta }}\in {\Re }^{3\times 3}$ are the constant inertia and stiffness matrices of vibration. On the other hand, other matrices $M_{\star },C_{\star }\in {\Re }^{3\times 3}$ are configuration-dependent nonlinear inertia and Coriolis matrices, with the off-diagonal terms representing the dynamic coupling.

Snapshots and profiles of kinetic, potential and total energies of undamped vibrations of this flexible beam simulated by our PMI-based framework with x=0 are shown in Figure 12 and Figure 13, respectively, where the integration interval is 1 kHz. There, we can see that the total energy is conserved throughout this fairly fast undamped vibration, clearly showing the ability of our PMI-based framework to simulate lossless, marginally stable or energy-conserving behavior. This undamped harmonic oscillation cannot be stably simulated by the SEI. It can however be done by the IEI, for which we must adopt some linearization as introduced in Duriez et al. (2006) and Courtecuisse et al. (2015) to render the IEI interactively fast. The results with this IEI simulation are shown in Figure 14, where it is clear that the IEI exhibits strictly stable oscillation with the total energy quickly decaying to zero. This lack of ability of the IEI to simulate lossless behavior is in fact due to its dissipativity (see Section 2.3), meaning that its advantage over our PMI framework is unclear, as the implementation complexities of PMI and IEI are at least comparable.

Fig. 12.

Snapshots of PMI-based simulation of undamped flexible beam harmonic oscillation.

Fig. 13.

Profiles of kinetic (KE), potential (PE), and total (KE+PE) energies of the PMI-based simulation of undamped flexible beam vibration in Figure 12: note the exact conservation of total energy during fast oscillations.

Fig. 14.

Profiles of kinetic (KE), potential (PE), and total (KE+PE) energies of the IEI-based simulation of undamped flexible beam vibration in Figure 12 with linearized inertia matrix: note the total energy quickly decays due to dissipativity of the IEI, although the same initial conditions as in Figure 13 are used.

For the haptic interaction, we connect the flexible beam to a square rigid box via a spring, which is in turn connected to an Omega 3 haptic device via virtual coupling (Lee et al., 2012); see Figure 15. The dynamics of the total system can then be written as

$$\begin{array}{c}\left[\begin{array}{cc}{M}_r\hfill & 0\hfill \\ 0\hfill & M(q)\hfill \end{array}\right]\left(\begin{array}{c}{ẍ}_r\hfill \\ \stackrel{··}{q}\hfill \end{array}\right)+\left[\begin{array}{cc}{B}_r\hfill & 0\hfill \\ 0\hfill & C(q)\hfill \end{array}\right]\left(\begin{array}{c}{ẋ}_r\hfill \\ \stackrel{.}{q}\hfill \end{array}\right)\\ +\left[\begin{array}{cc}{K}_r\hfill & -K_{\mathrm{rq}}\hfill \\ -K_{\mathrm{rq}}^T\hfill & K_q\hfill \end{array}\right]\left(\begin{array}{c}{x}_r\hfill \\ q\hfill \end{array}\right)=\left(\begin{array}{c}{K}_{\mathrm{vc}}{x}_{\mathrm{hd}}\hfill \\ f\hfill \end{array}\right)\end{array}$$

where $q:=[x;\delta ]\in {\Re }^6$ , and $x_r,x_{\mathrm{hd}}\in$ SE(2) are the maximal coordinates of the rigid box and the haptic device, $M_r=0.001I_{3\times 3}$ kg is the inertia of the rigid body, $K_{\mathrm{rq}}:=[K_{\mathrm{rq},x};0_{3\times 3}]\in {\Re }^{3\times 6}$ where $K_{\mathrm{rq},x}=10^5{I}_{3\times 3}$ is the spring between the beam end and the rigid box, $K_r:=K_{\mathrm{rq},x}+K_{\mathrm{vc}}\in {\Re }^{3\times 3}$ where $K_{\mathrm{vc}}=7500I_{3\times 3}$ is the virtual coupling spring, and $B_r$ is the stabilizing damping for the virtual spring. We also set the rotation of $x_{\mathrm{hd}}$ to be zero, since Omega 3 is the device of E(3).

Fig. 15.

Flexible beam (in generalized coordinates) is connected to rigid body (in maximal coordinates) via spring, which is in turn connected to haptic device via virtual coupling.

We then implement the virtual coupling according to Lee et al. (2012), which was derived to be applicable to PMI while also enforcing (hybrid-time) two-port passivity of the virtual coupling. For instance, with no virtual coupling damping, the force feedback to the haptic device is given by

$$\begin{array}{c}{f}_{\mathrm{hd}}(t):=-K_{\mathrm{vc}}(x_{\mathrm{hd},k}-x_{r,k}),t\in [t_k,t_{k+1})\hfill \end{array}$$

which, even if it contains $x_{\mathrm{hd},k}$ and $x_{r,k}$ instead of (continuous-time) $x_{\mathrm{hd}}(t)$ and (acausal) ${\widehat{x}}_{r,k}$ , can still be made passive by the device damping in the continuous-time domain and by the virtual damping $B_r$ in the discrete-time domain. The device damping of our adopted Omega 3 device turns out to be fairly small. This then requires us to choose a rather small force-scaling factor $\alpha =0.001$ . This small $\alpha$ is in fact an inherent restriction imposed by the device itself (i.e. z -width; Colgate and Brown, 1994). See Huang and Lee (2011) and Lee et al. (2012) for more details on the derivation and passivity of the virtual coupling for the PMI. Of course, we may choose larger (i.e. sharper) gain $K_{\mathrm{vc}}$ if we adopt techniques less conservative than the proportional–derivative-type virtual coupling (e.g. passive set-position modulation; Lee and Huang, 2010).

A profile of haptic feedback along the vertical z -direction of the haptic device during the interaction with the flexible beam is shown in Figure 16, where the user haptically interacts with the (lossless/marginally stable) flexible beam through the rigid box and virtual coupling, partly formulated in generalized coordinates and partly in maximal coordinates.

Fig. 16.

Haptic feedback from the flexible beam to the haptic device through virtual coupling connected to the rigid box.

5.3. UATD hand grasping

The last example is the adaptive power-grip of a virtual rigid ball by a multi-fingered UATD hand as shown in Figure 17. For the hand, we emulate the configuration of the iRobot-Harvard-Yale hand (Odhner et al., 2014). We model each finger to be constructed according to the design of Ozawa et al. (2014). For this case, it is then much more convenient to model each finger in generalized coordinates (i.e. joint angles), but model the hand’s palm and the virtual ball in maximal coordinates in SE(3).

Fig. 17.

The iRobot-Harvard-Yale hand (top left), the three-fingered virtual UATD hand used in the simulation (bottom left), and the tendon routing of each finger (right).

Each finger is of three DOF, but with only 2-DOF actuation via active tendons. See Figure 17 for their tendon routing (adopted from Ozawa et al., 2014). There are two active tendons (green and orange lines) and four passive tendons (other lines). The configuration of the UATD hand is then given by the nine joint angles $q\in {\Re }^9$ , which are actuated by the six actuators (i.e. motors) with their angles denoted by $\theta \in {\Re }^6$ . The total stiffness of the three fingers, stemming both from the active and the passive tendons, can then be computed more conveniently in the joint angles q (i.e. generalized coordinates) such that

$$K=J_a^T{K}_a{J}_a+J_p^T{K}_p{J}_p$$

where $K_a\in {\Re }^{6\times 6}$ , $K_p\in {\Re }^{12\times 12}$ are the active and passive tendon stiffnesses, and $J_a\in {\Re }^{6\times 9}$ , $J_p\in {\Re }^{12\times 9}$ are their respective Jacobians. Here, $K_a,K_p$ are typically assumed to be constant (i.e. linear tendons), whereas $J_a$ , $J_p$ are also constant, as they are defined by the pulley radii at each joint of the active and passive tendons, with $J_{a}q$ and $J_{p}q$ being their respective length changes. The dynamics of the UATD hand can then be written with respect to the generalized coordinates q:

$$\begin{array}{c}M(q)\stackrel{··}{q}+C(q,\stackrel{.}{q})\stackrel{.}{q}+\mathrm{Kq}=f+f_c\hfill \end{array}$$

(51)

similar to (4), where $M,C$ are the inertia and Coriolis matrices, Kq is the total tendon compliance, $f_c$ is the (external) grasping force, and f is the joint torque control input as given by

$$f=-J_a^T{K}_a{R}_a\theta =:-J_{\theta }\theta$$

where $R_a\in {\Re }^{6\times 6}$ is the matrix of the motor pulley radii (i.e. the tendon tension $K_a{R}_a\theta$ converted to the joint torque via $J_a^T$ ).

To control the adaptive under-actuated power-grip of this UATD hand, we use the Omega 3 haptic device as in Figure 15. Here, we only use one DOF of the haptic device (i.e. z -direction vertical motion) to command the power grip. For this, we first map $z\in [z_{\min },z_{\max }]\in \Re$ of the haptic device to the line interval $q^d\in \mathrm{line}(q_1^d,q_2^d)$ of the UATD hand configuration:

$$q^d(z):=q_1^d+\frac{z-z_{\min }}{z_{\max }-z_{\min }}(q_2^d-q_1^d)$$

where $q_1^d,q_2^d$ represent unfolded and folded configurations. Here, $q_2^d$ is not necessarily required to form the shape of the object to grasp, since the UATD hand can adaptively grip objects with different shapes via its tendon compliance and under-actuation. Then, using the quasi-static equation of (51), we can compute the following control input $\theta$ to attain $q^d(z)$ :

$$\theta =-{(J_{\theta }^T{J}_{\theta })}^{-1}{J}_{\theta }^{T}K{q}^d(z)$$

Of course, the UATD hand in general cannot exactly follow this $q^d(z)$ command, since its dynamics (51) is under-actuated; thus, the contact force $f_c$ can induce (unintended) deformation of the fingers. We provide this information to the user so that (s)he can adjust the adaptive grasping in the form of the following haptic device force feedback:

$$f_z(z,z_{\mathrm{proj}}(q)):=-k_{\mathrm{hd}}(z-z_{\mathrm{proj}}(q))$$

where $k_{\mathrm{hd}}>0$ is the feedback gain, z is the device’s current z -position, and $z_{\mathrm{proj}}(q)$ the mapping of the UATD hand configuration q to the z -direction of the haptic device:

$$z_{\mathrm{proj}}(q):=z_{\min }+(z_{\max }-z_{\min })\frac{\parallel q-q_1^d\parallel }{\parallel q_2^d-q_1^d\parallel }$$

Snapshots of this adaptive power-grip of the multi-fingered UATD hand are shown in Figure 18, and the contact force feedback is shown in Figure 19. For this, we use $K_a=100I_{6\times 6}$ N/m and $K_p=10I_{12\times 12}$ N/m with the mass of each link of the finger being 0.001 kg. From Figure 18 and Figure 19, we can then see that our PMI-based framework can stably simulate this power-gripping of the UATD hand, through the PMI and the multi-point PMI-LCP formulated both for the fingers (i.e. in generalized coordinates: Section 4) and for the virtual ball and palm (i.e. in maximal coordinates: Section 3). Note also that the fingers are fairly light/stiff systems, which are difficult to stably simulate by using other techniques (e.g. SEI).

Fig. 18.

Snapshots of adaptive power-grasping of virtual rigid ball by the three-fingered UATD hand, with each finger having tendon-routing as depicted in Figure 17, with the three fingers simultaneously controlled by one DOF of the haptic device.

Fig. 19.

Haptic feedback of the 1-DOF (i.e. z-directional) contact force of the adaptive power-grasping of the UATD hand in Figure 18.

6. Summary and future research

In this paper, we develop a novel framework for haptic rendering and interactive simulation which is based on the PMI, which is derived by directly enforcing the passivity property in the discrete-time domain. We formulate this PMI-based framework for mechanical systems, both in maximal coordinates and generalized coordinates; we also devise a PMI-LCP formulation to fully incorporate multi-point Coulomb friction contact modeling into these two coordinate settings. Thanks to its passivity property, our PMI-based framework allows us:

1. to tune the virtual objects freely by exploring a wide range of parameters in a modular manner;

2. to stably simulate very light/stiff objects (e.g. for transparency) even with variable/slow integration steps (e.g. due to intermittent multi-point contact); and

3. to emulate such marginally stable, lossless or energy-conserving behavior as the harmonic oscillation of flexible objects or elastic contact between rigid objects.

Furthermore, our PMI-based framework provides implementation flexibility and convenience, by allowing for both maximal and generalized coordinates, and also their combination, which are imperative when simulating some types of mechanical systems (e.g. a multi-fingered UATD hand). These advantages are experimentally demonstrated in this paper by using three illustrative examples: haptic operation of a peg-in-hole task; haptic interaction with a flexible beam; and adaptive power-grasping by a multi-fingered UATD hand.

Some possible future research directions include the following: (1) Extension of the PMI-based framework to large-scale deformable objects, including cutting-and-joining and fluidic simulation; (2) Combining the constraint-stabilizing potential and PMI-LCP formulation while respecting non-iterativeness and passivity of the PMI; (3) Geometric properties of the coordinate transformation (36) and its application to other problems; (4) Application of the PMI-based rendering for multi-user virtual reality and haptic interaction over the Internet; (5) Development of an open-source simulation platform based on our PMI-based simulation framework.