System and method for executing a simulation of a constrained multi-body system

Abstract

Methods of and systems for executing a simulation of a constrained multi-body system. The method comprises, using a physics engine, simulating the constrained multi-body system, wherein: the constrained multi-body system comprises articulated constraints, the articulated constraints are associated with a geometric stiffness matrix; the geometric stiffness matrix defining a geometric stiffness; a diagonal approximation of the geometric stiffness matrix is generated; and the diagonal approximation is used as part of a stability analysis in which damping is automatically adjusted so that the damping stabilizes the simulation of the constrained multi-body system.

Claims

1. A computer-implemented method for executing a simulation of a constrained multi-body system, the method comprising: simulating the constrained multi-body system by using a physics engine, the constrained multi-body system comprising articulated constraints; associating the articulated constraints with a geometric stiffness matrix; the geometric stiffness matrix defining a geometric stiffness, wherein the geometric stiffness defines how a direction of forces in the constrained multi-body system changes due to instantaneous motion of one or more links associated with the constrained multi-body system; and using the geometric stiffness matrix as part of a stability analysis in order to determine damping values, the damping values being calculated at each time step of the simulation, the damping values being determined so as to achieve stability of the simulation of the constrained multi-body system.

2. The method of claim 1, wherein determining the damping values is based on a damping coefficient.

3. The method of claim 2, wherein the damping coefficient is automatically adjusted so as to be large enough to stabilize transverse oscillations of the simulation of the constrained multi-body system.

4. The method of claim 1, wherein the simulation of the constrained multi-body system comprises a dynamical system wherein the geometric stiffness matrix is reinterpreted as a spring and the generated forces are used in combination with a stability criterion for stable simulation.

5. The method of claim 4, wherein spring forces define the damping.

6. The method of claim 4, wherein the stability criterion comprises a threshold based on a behavior of a simple 1D undamped oscillator.

7. The method of claim 6, wherein the geometric stiffness matrix serves as the stability criterion as to whether the dynamical system becomes unstable.

8. The method of claim 7, wherein the geometric stiffness matrix is a {tilde over (K)} matrix, wherein the {tilde over (K)} matrix is defined as follows: $\tilde{K}=\frac{J^T}{x}$ wherein J is a matrix of constraint Jacobians, are constraint forces and x are positions and orientations of bodies in the simulation of the constrained multi-body system.

9. The method of claim 1, wherein the instantaneous motion comprises at least one of rotation of the one or more links associated with the constrained multi-body system, translation of the one or more links associated with the constrained multi-body system and relative motion of the one or more links associated with the constrained multi-body system.

10. The method of claim 1, wherein a topology of a mechanical structure of the constrained multi-body system remains unmodified during the simulating of the constrained multi-body system.

11. The method of claim 1, wherein the damping values comprise parameters which define an interpretation of the geometric stiffness used in the simulation of the constrained multi-body system.

12. The method of claim 11, wherein at least one of the parameters define physical behaviour of a transient physical spring that is able to generate forces in directions transverse to transverse oscillations directions.

13. The method of claim 12, wherein data encoding parameters of the transient spring models a tensor encoding variations in constraint force directions.

14. The method of claim 1, wherein the simulating of the constrained multi-body system is executed in real-time.

15. The method of claim 1, wherein stabilizing the simulation of the constrained multi-body system comprises stabilizing the transverse oscillations of the simulation of the constrained multi-body system.

16. The method of claim 1, wherein the geometric stiffness matrix is generated by using a low rank update algorithm, the using of the low rank update algorithm comprising constructing a lead matrix of dynamical equations of motion, the lead matrix comprising the geometric stiffness matrix.

17. The method of claim 1, wherein the geometric stiffness matrix models the articulated constraints.

18. The method of claim 1, wherein using the geometric stiffness matrix comprises generating a diagonal approximation of the geometric stiffness matrix.

19. The method of claim 18, wherein the diagonal approximation of the geometric stiffness is calculated in accordance with the following equation: $b=\frac{h^{2}k-4m}{2h}$ wherein k is a stiffness, b is a damping and m is a mass, h is a given time step and a is a positive scalar value that allows control over how close dynamical system may be to the boundary before damping may be required.

20. The method of claim 18, wherein the diagonal approximation of the geometric stiffness matrix preserves mechanical work properties associated with the constrained multi-body system.

21. A computer-implemented system for executing a simulation of a constrained multi-body system, the system comprising: a processor; a non-transitory computer-readable medium, the non-transitory computer-readable medium comprising control logic which, upon execution by the processor, causes: simulating the constrained multi-body system by using a physics engine, the constrained multi-body system comprising articulated constraints; associating the articulated constraints with a geometric stiffness matrix; the geometric stiffness matrix defining a geometric stiffness, wherein the geometric stiffness defines how a direction of forces in the constrained multi-body system changes due to instantaneous motion of one or more links associated with the constrained multi-body system; and using the geometric stiffness matrix as part of a stability analysis in order to determine damping values, the damping values being calculated at each time step of the simulation, the damping values being determined so as to achieve stability of the simulation of the constrained multi-body system.

22. A non-transitory computer-readable medium, the non-transitory computer-readable medium comprising control logic which, upon execution by a processor, causes: simulating the constrained multi-body system by using a physics engine, the constrained multi-body system comprising articulated constraints; associating the articulated constraints with a geometric stiffness matrix; the geometric stiffness matrix defining a geometric stiffness, wherein the geometric stiffness defines how a direction of forces in the constrained multi-body system changes due to instantaneous motion of one or more links associated with the constrained multi-body system; and using the geometric stiffness matrix as part of a stability analysis in order to determine damping values, the damping values being calculated at each time step of the simulation, the damping values being determined so as to achieve stability of the simulation of the constrained multi-body system.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) For a better understanding of the present technology, as well as other aspects and further features thereof, reference is made to the following description which is to be used in conjunction with the accompanying drawings, where:

(2) FIG. 1 is a diagram of a computing environment in accordance with an embodiment of the present technology;

(3) FIG. 2 is a diagram of rigid bodies in accordance with embodiments of the present technology;

(4) FIG. 3 is a representation of a simulation and a graph in accordance with embodiments of the present technology;

(5) FIG. 4 is a representation of graphs in accordance with embodiments of the present technology;

(6) FIG. 5 is a representation of a graph in accordance with embodiments of the present technology;

(7) FIG. 6 is a representation of a simulation and a graph in accordance with embodiments of the present technology;

(8) FIG. 7 is a representation of a graph in accordance with embodiments of the present technology;

(9) FIG. 8 is a representation of a simulation and a graph in accordance with embodiments of the present technology;

(10) FIG. 9 is a representation of a simulation in accordance with embodiments of the present technology;

(11) FIGS. 10 and 11 are a representation of a simulation in accordance with embodiments of the present technology; and

(12) FIG. 12 is a diagram illustrating a flowchart illustrating a computer-implemented method implementing embodiments of the present technology.

(13) It should also be noted that, unless otherwise explicitly specified herein, the drawings are not to scale.

DETAILED DESCRIPTION

(14) The examples and conditional language recited herein are principally intended to aid the reader in understanding the principles of the present technology and not to limit its scope to such specifically recited examples and conditions. It will be appreciated that those skilled in the art may devise various arrangements which, although not explicitly described or shown herein, nonetheless embody the principles of the present technology and are included within its spirit and scope.

(15) Furthermore, as an aid to understanding, the following description may describe relatively simplified implementations of the present technology. As persons skilled in the art would understand, various implementations of the present technology may be of a greater complexity.

(16) In some cases, what are believed to be helpful examples of modifications to the present technology may also be set forth. This is done merely as an aid to understanding, and, again, not to define the scope or set forth the bounds of the present technology. These modifications are not an exhaustive list, and a person skilled in the art may make other modifications while nonetheless remaining within the scope of the present technology. Further, where no examples of modifications have been set forth, it should not be interpreted that no modifications are possible and/or that what is described is the sole manner of implementing that element of the present technology.

(17) Moreover, all statements herein reciting principles, aspects, and implementations of the present technology, as well as specific examples thereof, are intended to encompass both structural and functional equivalents thereof, whether they are currently known or developed in the future. Thus, for example, it will be appreciated by those skilled in the art that any block diagrams herein represent conceptual views of illustrative circuitry embodying the principles of the present technology. Similarly, it will be appreciated that any flowcharts, flow diagrams, state transition diagrams, pseudo-code, and the like represent various processes which may be substantially represented in computer-readable media and so executed by a computer or processor, whether or not such computer or processor is explicitly shown.

(18) The functions of the various elements shown in the figures, including any functional block labeled as a processor, may be provided through the use of dedicated hardware as well as hardware capable of executing software in association with appropriate software. When provided by a processor, the functions may be provided by a single dedicated processor, by a single shared processor, or by a plurality of individual processors, some of which may be shared. In some embodiments of the present technology, the processor may be a general purpose processor, such as a central processing unit (CPU) or a processor dedicated to a specific purpose, such as a digital signal processor (DSP). Moreover, explicit use of the term a processor, a simulation application should not be construed to refer exclusively to hardware capable of executing software, and may implicitly include, without limitation, application specific integrated circuit (ASIC), field programmable gate array (FPGA), read-only memory (ROM) for storing software, random access memory (RAM), and non-volatile storage. Other hardware, conventional and/or custom, may also be included.

(19) Software modules, or simply modules which are implied to be software, may be represented herein as any combination of flowchart elements or other elements indicating performance of process steps and/or textual description. Such modules may be executed by hardware that is expressly or implicitly shown. Moreover, it should be understood that module may include for example, but without being limitative, computer program logic, computer program instructions, software, stack, firmware, hardware circuitry or a combination thereof which provides the required capabilities.

(20) With these fundamentals in place, we will now consider some non-limiting examples to illustrate various implementations of aspects of the present technology.

(21) FIG. 1 illustrates a diagram of a computing environment 100 in accordance with an embodiment of the present technology is shown. In some embodiments, the computing environment 100 may be implemented by any of a conventional personal computer, a computer dedicated to operating simulation and/or an electronic device (such as, but not limited to, a mobile device, a video game system, a tablet device, etc.) and/or any combination thereof appropriate to the relevant task at hand. In some embodiments, the computing environment 100 comprises various hardware components including one or more single or multi-core processors collectively represented by a processor 110 a solid-state drive 120, a random access memory 130 and an input/output interface 150. The computing environment 100 may be a computer specifically designed for operating simulation engines and/or physics engines. In some alternative embodiments, the computing environment 100 may be a generic computer system.

(22) In some embodiments, the computing environment 100 may also be a sub-system of one of the above-listed systems. In some other embodiments, the computing environment 100 may be an off the shelf generic computer system. In some embodiments, the computing environment 100 may also be distributed amongst multiple systems. The computing environment 100 may also be specifically dedicated to the implementation of the present technology. As a person in the art of the present technology may appreciate, multiple variations as to how the computing environment 100 is implemented may be envisioned without departing from the scope of the present technology.

(23) Communication between the various components of the computing environment 100 may be enabled by one or more internal and/or external buses 160 (e.g. a PCI bus, universal serial bus, IEEE 1394 Firewire bus, SCSI bus, Serial-ATA bus, ARINC bus, etc.), to which the various hardware components are electronically coupled.

(24) The input/output interface 150 may allow enabling networking capabilities such as wire or wireless access. As an example, the input/output interface 150 may comprise a networking interface such as, but not limited to, a network port, a network socket, a network interface controller and the like. Multiple examples of how the networking interface may be implemented will become apparent to the person skilled in the art of the present technology. For example, but without being limitative, the networking interface may implement specific physical layer and data link layer standard such as Ethernet, Fibre Channel, Wi-Fi or Token Ring. The specific physical layer and the data link layer may provide a base for a full network protocol stack, allowing communication among small groups of computers on the same local area network (LAN) and large-scale network communications through routable protocols, such as Internet Protocol (IP).

(25) According to implementations of the present technology, the solid-state drive 120 stores program instructions suitable for being loaded into the random access memory 130 and executed by the processor 110 for executing a simulation of a constrained multi-body system. For example, the program instructions may be part of a library or an application.

(26) Notation

(27) In FIG. 2, a diagram 200 illustrates coordinate frames and offsets used to specify configuration of a joint 202. A body 204 (also referred to as the body i) and a body 206 (also referred to as the body j) are coupled by the joint 202 with attachment offsets s.sub.i, s.sub.j, coordinate frames {u.sub.i, v.sub.i, n.sub.i} and {u.sub.j, v.sub.j, n.sub.j} and displacement vector d.sub.ij. The attachment basis and offsets are fixed to the body 204 and the body 206. The offset relative to the center of mass of body i is specified by s.sub.i, which is represented in a global coordinate frame but specified by the local vector s.sub.i such that s.sub.i=R.sub.i s.sub.i. Similarly, orthogonal basis {u.sub.i, v.sub.i, n.sub.i} are vectors in the global coordinate frame, but specified by a local basis {u.sub.i, v.sub.i, n.sub.i} and transformed by R.sub.i. The dependence on a body orientation is relevant for deriving geometric stiffness matrices, such as the geometric stiffness matrices detailed below in connection with the section entitled Constraint library. Unless otherwise stated vector quantities are represented in a global coordinate frame. Table 1 illustrated below provides an overview of variables and notation used throughout the present document.

(28) TABLE-US-00001 TABLE 1 Variables and notation used throughout this document. p.sub.i position of body i .sub.i orientation of body i {dot over (p)}.sub.i linear velocity of body i .sub.i angular velocity of body i x.sub.i spatial configuration of body i, such that x.sub.i = (p.sub.i.sup.T .sub.i.sup.T).sup.T x configuration of all bodies as (x.sub.0.sup.T x.sub.1.sup.T . . . x.sub.n.sup.T).sup.T v velocity of all bodies as ({dot over (p)}.sub.0.sup.T .sub.0.sup.T . . . {dot over (p)}.sub.n.sup.T .sub.n.sup.T).sup.T R.sub.i rotation matrix corresponding to the orientation .sub.i s.sub.i joint attachment offset in local frame s.sub.i joint attachment offset in the global frame, s.sub.i = R.sub.is.sub.i a position level constraint equation J $a\mathrm{constraint}\mathrm{Jacobian}\mathrm{matrix},\mathrm{such}\mathrm{that}J=\frac{d}{\mathrm{dx}}$ J Jacobian matrix for all constraints {tilde over (K)} geometric stiffness matrix {circumflex over ()} skew-symmetric cross product matrix of a vector .sub.+ an implicit term or quantity

(29) Geometric Stiffness

(30) In one aspect of the present technology, geometric stiffness may be defined as a tensor encoding variations in constraint force directions. In some embodiments, the tensor may have the following form:

(31) $\begin{array}{cc}\tilde{K}=\frac{J^T}{x}.& \left(1\right)\end{array}$

(32) In this equation, J is a matrix of constraint Jacobians, are constraint forces and x are positions and orientations of bodies in a simulation. Tournier et al. (TOURNIER M., NESME M., GILLES B., FAURE F.: Stable constrained dynamics. ACM Transactions on Graphics (TOG) 34, 4 (2015), 132.), the entirety of which is hereby incorporated by reference in jurisdictions where such incorporation is permissible, introduces the term geometric stiffness as an implicit stiffness in the velocity-level discretization of the constrained dynamical equations:

(33) $\begin{array}{cc}\left(\begin{array}{cc}\tilde{M}& -J^T\\ J& \frac{1}{h^2}C\end{array}\right)\left(\begin{array}{c}{v}_{+}\\ h{}_{+}\end{array}\right)=\left(\begin{array}{c}p+\mathrm{hf}\\ -\frac{1}{h}\end{array}\right)& \left(2\right)\end{array}$

(34) In this equation, {tilde over (M)}=Mh.sup.2{tilde over (K)}, p=Mv is a current momentum, f contains external forces, diagonal regularization matrix C makes constraints compliant and h is the time step size. Forming a Schur complement of the upper left block of the equation gives the following reduced system:

(35) $\begin{array}{cc}\left(\frac{1}{h^2}C+J{\tilde{M}}^{-1}{J}^T\right)h+=-\frac{1}{h}-J{\tilde{M}}^{-1}\left(p+\mathrm{hf}\right)& \left(3\right)\end{array}$

(36) The reduced system above may be similar to the one used by some open source and/or commercial rigid body physics engines (such as, but not limited to, Bullet, Open Dynamics Engine, Havok, Simmechanics from Matlab and Vortex from CM Labs). The reduced system above typically requires solving a mixed linear complementarity problem (MLCP) as the system may include both bilateral and unilateral constraints. As it may be appreciated by the person skilled in the art of the present technology, inclusion of geometric stiffness in equations (2) and/or (3) may lead to numerical difficulties upon solving linear systems. In some instances, such difficulties may relate to symmetry, positive definiteness and/or efficiency.

(37) With respect to difficulties relating to symmetry and positive definiteness, in some instances, solving the dense linear system (3) may involve a Cholesky factorization. This may require {tilde over (M)} to be symmetric and positive definite. This may also be a similar requirement for the system (2). However, inclusion of {tilde over (K)} may result in a non-symmetric matrix. Tournier el al. suggests that a symmetric approximation of the matrix may be formed with little effect on stabilizing properties. However, the resulting matrix may still fail to meet positive definite requirement of the numerical method. Alternative approaches comprise Teran et al. (TERAN J., SIFAKIS E., IRVING G., FEDKIW R.: Robust quasistatic finite elements and flesh simulation. In Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation (2005), ACM, pp. 181-190.).), the entirety of which is hereby incorporated by reference in jurisdictions where such incorporation is permissible, which describes a method for enforcing positive definiteness by clamping eigenvalues of a diagonalized elastic deformation gradient to be non-negative. This method may lead to robust simulations but can be costly and therefore may make the method not well suited for real-time and/or performance-focused applications.

(38) With respect to difficulties relating to efficiency, solving the reduced system of equation (3) may require forming an inverse of {tilde over (M)}. For a block diagonal mass matrix, which may be typically the case, this may be done efficiently by inverting linear mass and performing a fast 33 inversion of an inertia of each body. However, with inclusion of geometric stiffness, inverting the full mass matrix may no longer be straightforward due to nonzero terms that may appear outside the block diagonal.

(39) To, at least partially, overcome some of the above-discussed difficulties, the present technology provides alternate formulations of geometric stiffness and inclusion of the alternate formulations of the geometric stiffness in multibody equations. As further discussed in the paragraphs below, the inverse of {tilde over (M)} may be obtained more efficiently by using low rank updates. In addition, the present technology also provides using diagonalization of the geometrix stiffness matrix as part of a stability analysis. This aspect with also be discussing in the paragraphs below.

(40) Inverse by Low Rank Updates

(41) Some aspects of the present technology relate to a method of constructing an inverse of {tilde over (M)} using low rank updates. In some circumstances, an unmodified matrix may be relatively efficiently inverted due to a nearly diagonal form of M. However, inverting an augmented mass matrix may require more computational effort (e.g., a lower upper (LU) decomposition) due to off-diagonal elements in a geometrical stiffness matrix.

(42) It should be appreciated that in many instances, geometric stiffness matrix associated with articulated rigid body constraints are sparse and low rank. In some embodiments of the present technology, these properties may be exploited to build an augmented mass matrix more efficiently. In some embodiments, an approach to build the augmented mass matrix may be based on the Sherman-Morrison formula as described in Sherman et al. (SHERMAN J., MORRISON W. J.: Adjustment of an inverse matrix corresponding to a change in one element of a given matrix. The Annals of Mathematical Statistics 21, 1 (1950), 124-127.), the entirety of which is hereby incorporated by reference in jurisdictions where such incorporation is permissible.

Example 1: Decomposition for a Simple Joint

(43) This first example aims at presenting the decomposition associated with a simple articulated joint, in this example, a ball-and-socket. In this example, a geometric stiffness matrix associated with the simple articulated joint is as follows:

(44) $\begin{array}{cc}{\tilde{K}}_{\mathrm{bs}}=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& \hat{}{\hat{s}}_i& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& -\hat{}{\hat{s}}_j\end{array}\right){}^{1212}& \left(4\right)\end{array}$

(45) In the geometric stiffness matrix (4), .sup.3 are the Lagrange multipliers enforcing the constraint. As it may be appreciated by the person skilled in the art of the present technology, the geometric stiffness matrix (4) may be qualified as being sparse and associated with a simple structure. In some embodiments of the present technology, a rank decomposition of {tilde over (K)}.sub.bs may be found by examining one of the non-zero blocks of (4) and identifying that s is in a right null space of the matrix {circumflex over ()}.

(46) Starting with a unit vector

(47) $\stackrel{.fwdarw.}{s}=\frac{s}{||s||}{}^3,$
an orthonormal basis may be constructed where:
{right arrow over (s)}{right arrow over (t)},{right arrow over (r)}{right arrow over (t)},{right arrow over (r)}{right arrow over (s)},
|{right arrow over (s)}|=|{right arrow over (t)}|=|{right arrow over (r)}|=1.

(48) A column space of {circumflex over ()} spans a subspace formed by vectors {right arrow over (r)} and {right arrow over (t)}. A vector of Lagrange multipliers may be projected onto this basis, giving the following coefficients:
.sub.r=.sup.T{right arrow over (r)}
.sub.s=.sup.T{right arrow over (s)}
.sub.t=.sup.T{right arrow over (t)}.

(49) As it may be appreciated, {right arrow over (r)}, {right arrow over (s)}, and {right arrow over (t)} are in a row space of the matrix {circumflex over ()}. Its decomposition is therefore a rank-2 matrix which is as follows:
{circumflex over ()}=s(.sub.r{right arrow over (s)}{right arrow over (r)}.sup.T.sub.s{right arrow over (r)}{right arrow over (r)}.sup.T+.sub.t{right arrow over (s)}{right arrow over (t)}.sup.T.sub.s{right arrow over (t)}{right arrow over (t)}.sup.T)(5)

(50) This decomposition (5) may be applied to both non-zero 33 blocks in the geometric stiffness matrix (4) and used to reconstruct the full matrix {tilde over (K)}.sub.bs. Knowing a rank decomposition of the matrix {tilde over (K)}.sub.bs and being able to compute an inverse mass matrix M.sup.1, if may become possible to compute an inverse of {tilde over (M)} by applying the Sherman-Morrison formula.

(51) The Sherman-Morrison Formula

(52) Given a rank-1 matrix ab.sup.T the Sherman-Morrison formula may be used to update the inverse mass matrix M.sup.1 by:

(53) ${\left(M+{\mathrm{ab}}^T\right)}^{-1}=M^{-1}-\frac{M^{-1}{\mathrm{ab}}^T{M}^{-1}}{1+b^T{M}^{-1}a}.$

(54) An augmented mass matrix {tilde over (M)}.sup.1 may therefore be obtained by sequentially applying rank-1 updates using a geometric stiffness decomposition of all joints in the system.

(55) As an example, we may consider an update to the block of body i using .sub.r.sub.i{right arrow over (s)}.sub.i{right arrow over (r)}.sub.i.sup.T, which is the first term in (5). We may define vectors , b.sup.6n with dimensionality matching an n body system. These vectors are zero except entries:
a6i+3 . . . 6i+5=.sub.r.sub.i{right arrow over (s)}.sub.i
b6i+3 . . . 6i+5={right arrow over (r)}.sub.i.

(56) There vectors are the 3D basis vectors padded with zeros according to block's location in the system matrix. A low rank update of the inverse mass matrix M.sup.6n6n may then be performed by:

(57) $\begin{array}{cc}{M}^{-1}{M}^{-1}+h^2||s||\frac{M^{-1}\stackrel{\_}{a}{\stackrel{\_}{b}}^T{M}^{-1}}{1+{\stackrel{\_}{b}}^T{M}^{-1}\stackrel{\_}{a}}.& \left(6\right)\end{array}$

(58) The scaling factors h.sup.2 and s are based on the inclusion of the geometric stiffness in the dynamical equations and length of the attachment offset s. Similar updates may be performed using the rank-1 matrices {right arrow over (s)}{right arrow over (t)}.sup.T, {right arrow over (r)}{right arrow over (r)}.sup.T, and {right arrow over (t)}{right arrow over (t)}.sup.T and for both bodies. In an embodiment of the present technology, the matrix-vector products of (6) may be computed by exploiting sparsity of the vectors for efficiency purposes. In some embodiments of the present technology, a method to incrementally construct an inverse augmented mass matrix may be similar for other types of joints. Low rank decompositions for joints used in many articulated rigid body simulators may be found in the section Examples of low rank decompositions for joints used in articulated rigid body simulators of the present document.

(59) It should be noted that the linear systems of (2) and (3) may not be positive all the time.

(60) Geometric Stiffness Inspired Damping

(61) The present section details a method for diagonalizing {tilde over (K)} so as to result in a positive definite matrix for the linear system. The diagonalized approximation may then be used as part of a stability analysis where damping may be introduced. In some embodiments, the damping may only be introduced when there is an indication that forces related to the geometric stiffness may cause instability.

(62) In some embodiments of the present technology, {tilde over (K)} may not represent a material property of the physical system but it may nonetheless be convenient to interpret a geometric stiffness as a physical spring that is able to generate forces in directions transverse to constraints. In some embodiments, the interpretation may facilitate analysis of the geometry stiffness as part of a mass-spring system of undamped harmonic oscillators.

(63) In order to illustrate an embodiment of the present technology, an example of a simple harmonic oscillator system will be detailed. Stability requirements for simulating the simple harmonic oscillator system is considered with a semi-implicit integration scheme which may also be referred to as symplectic Euler. The stability requirements may then be applied to an integrator with less strict stability requirements, such as a single step implicit integrator that is commonly used by interactive computer graphics applications. A stability criterion for stable simulation is derived and is extended to a dynamical system where spring forces are generated by the geometric stiffness, with the {tilde over (K)} matrix serving as an indicator of when the system may become instable. Thus, damping may be introduced to stabilize the system using an adaptive method (i.e., an adaptive damping scheme) which estimates a damping coefficient that may be large enough to stabilize the system.

(64) Stability Analysis of Harmonic Oscillator

(65) As an example, we consider a behavior of a 1D Hookean mass-spring-damper attached to a particle. Accelerations generated by such a dynamical system may be as follows:
{umlaut over (x)}=(kxb{dot over (x)})/m(7)

(66) Where x is a displacement from rest length, {dot over (x)} is a particle velocity, k is a stiffness, b is a damping and m is a mass of the particle. Simulating the dynamical system with a semi-implicit integrator and time step h gives the following velocity and position update:
{dot over (x)}{dot over (x)}+h{umlaut over (x)},
xx+h{dot over (x)}.

(67) Substituting (7) and representing updates in matrix form gives the following space equation:

(68) $\begin{array}{cc}\left(\begin{array}{c}x\\ \stackrel{.}{x}\end{array}\right)\underset{\underset{P}{}}{\left(\begin{array}{cc}1-h^2\frac{k}{m}& h-h\frac{b}{m}\\ -h\frac{k}{m}& 1-h\frac{b}{m}\end{array}\right)}\left(\begin{array}{c}x\\ \stackrel{.}{x}\end{array}\right)& \left(8\right)\end{array}$

(69) With the eigenvalues of matrix P given by:

(70) 0$-\frac{h^{2}k+\mathrm{bh}-2m\sqrt{h^4{k}^2+2{\mathrm{bh}}^{3}k+b^2{h}^2-4h^2\mathrm{km}}}{2m}$

(71) If the magnitude of either eigenvalue is greater than 1, then the simulation may become unstable.

(72) As it may be appreciated from a review of (8), multiple parameters may affect simulation stability. Mass, stiffness and time step may be typically fixed due to application or model requirements. Therefore, in some embodiments of the present technology, only damping may be used to stabilize the dynamical system. In some embodiments of the present technology, non-constitutive damping may be used to stabilize the simulations. In some approaches, non-constitutive damping may be tuned by manually tuning damping coefficients. This may be cumbersome since values must be selected that minimally affect the dynamics, yet leave the simulation stable for a wide range of configurations. To address, at least some of the drawbacks of the manual approach, some embodiments of the present technology rely on an automatic method for computing damping based on a stability region of a semi-implicit integrator.

(73) In some embodiments of the present technology, a single step backward Euler method with a linear approximation of the constraint Jacobian and mass matrix is used. The backward Euler method may bring stability to the simulation with a region of stability which may contain that of the semi-implicit integrator used in the analysis set forth above. Semi-implicit Euler may exhibit absolute stability when the eigenvalues of a phase space linear system lie within a negative real half of a complex plane. This is the region of stability for the integration method whereas backward Euler may have absolute stability for all but a small region in a positive real half of the complex plane. In some embodiments of the present technology, the present technology proposes a stability criterion based on an integrator with a smaller region of stability and then apply it to an implicit integrator. In some embodiments, the stability criterion may come in the form of a threshold based on a behavior of a simple ID undamped oscillator.

(74) Stability Threshold

(75) In some embodiments, a stability criterion (also referred to as a stability threshold) may be determined by setting b=0 in eigenvalue equations, giving an undamped oscillator. This may result in the following stability inequality:

(76) $\begin{array}{cc}\frac{h^{2}k-\sqrt{h^4{k}^2-4h^2\mathrm{km}}}{m}4,& \left(9\right)\end{array}$

(77) The stability inequality (9) may be true of eigenvalue magnitudes are less than 1.

(78) Alpha Parameter

(79) In some embodiments of the present technology, the stability inequality (9) may be augmented with a parameter which is a positive scalar value that allows control over how close dynamical system may be to the boundary before damping may be required. This gives the following inequality:

(80) $\begin{array}{cc}\frac{h^{2}k}{m}4.& \left(10\right)\end{array}$

(81) In some embodiments, if the inequality (10) is violated, then the damping required may be found by solving for b at the boundary. This may be computed as:

(82) $\begin{array}{cc}b=\frac{h^{2}k-4m}{2h}.& \left(11\right)\end{array}$

(83) For a spring system with constant mass and stiffness, the value of b may also be constant. However, if the spring parameters are estimated from the geometric stiffness, the damping values may change at each time step. The stability threshold may also vary with the constraint forces. In some embodiments, for the purpose of stabilizing a general class of physics simulation, a stability threshold may be evaluated at each time step and damping values may be computed according to the equation (11) thereby providing an adaptive damping scheme.

(84) In some embodiments, the above analysis may be extended to the {tilde over (K)} matrix. In some embodiments, the geometric stiffness is associated with the {tilde over (K)} matrix. Therefore, in some embodiments, a geometric stiffness matrix may be referred to as the {tilde over (K)} matrix. In some instances, modal analysis may require finding eigenvalues of a matrix M.sup.1{tilde over (K)} which may be computationally costly. As a result, in some embodiments, the present technology proposes an approximation of the matrix thereby rendering the stability analysis less costly computationally-wise. If the {tilde over (K)} matrix is diagonal, then the damping scheme for stabilizing a simple harmonic oscillator may be applicable.

(85) Diagonalizing the {tilde over (K)} Matrix

(86) The present section provides details on how a diagonal approximation simplifying a stability analysis may be computed, in accordance with embodiments of the present technology. One objective of the approach may be to compute a simplified matrix for use as a stability criterion. In some embodiments, the geometric stiffness may be interpreted as a transient spring. In some embodiments, mechanical work done by the {tilde over (K)} matrix may be considered in building a diagonal approximation K.sub.d. An instantaneous transfer of mechanical energy occurring over a single time step, or also referred to as instantaneous work, done by the geometric stiffness may be 0.5 v.sup.Tv.

(87) In some embodiments, a diagonal approximation may represent a spring that tends to do at least a same amount of mechanical work, by providing a stiffer system in directions of any typical v. An approximation that tends to upper bound the work may be computed by assigning each diagonal element K.sub.di,i norm of a corresponding column vector in {tilde over (K)}, such that:
K.sub.di,i=|{tilde over (K)}.sub.i|.(12)

(88) Once a diagonal approximation of the geometric stiffness matrix is generated, a stability analysis may be conducted more easily. In some embodiments, in an inertial coordinate frame, it may become a matter of assigning mass and stiffness values and apply equation (11) to compute entry B.sub.i,i of the diagonal damping matrix, the mass and stiffness values being as follows:
m=M.sub.i,i
k=K.sub.di,i

(89) In some embodiments, this may be done for all degree-of-freedoms (DOFs), although, in some embodiments, damping coefficients corresponding to translational DOFs may be (but not necessarily) sometimes discarded. An augmented mass matrix may then be rebuilt as the following equation:
{tilde over (M)}=M+hB.sub.;(13)

(90) Equation (13) may be positive definite since M may be positive definite and B may contain only non-negative values along the diagonal. In this embodiment, B may be seen as replacing the geometric stiffness matrix. As the damping coefficients may augment the mass and inertia of bodies in the system, it may also be referred to the above approach as an inertial damping formulation.

(91) Projection to Constraint Space

(92) In some alternative embodiments, damping may be applied to constraint rows of a linear system. Such approach may stabilize in directions that may be aligned with constraint axes. In some cases, simulations may appear more plausible when damping is applied to constraint equations. This may be achieved by mapping the mass and diagonalized geometric stiffness matrix to a constraint space and performing a stability analysis there. The mass M and the diagonal approximation K.sub.d may be mapped onto each constraint row by:
m=(JM.sup.1j.sup.T).sup.1
k=(JK.sub.d.sup.1J.sup.T).sup.1,

(93) Where J is the Jacobian for a single constraint equation and m and k are effective mass and stiffness, respectively, for the constraint. Coupling of constraints through the mass and stiffness matrix may be ignored. In some embodiments, rather than form an effective mass and stiffness using a full constraint matrix, J may simply be a row vector. This may result in a single scalar value for m and k. Stability analysis may be performed using m and k by applying the equation (11) to compute a damping coefficient, bfor the constraint. However, a constraint damping term may not be considered by the dynamical equations (2). In such circumstances, the formulation may be modified to support constraint damping.

(94) Constraint Space Damping

(95) In some embodiments, constraint equations may be reorganized as a force equation for a spring-damper system to define the following equation:
.sub.+=C.sup.1.sub.+B.sub.+(14)

(96) Where B may be a diagonal matrix of non-negative damping coefficients which may be assembled from the b computed for each constraint row. Substituting Jv.sub.+=.sub.+ and letting .sub.+=+h.sub.+, the constraint equations may become:

(97) $\begin{array}{cc}{\mathrm{Jv}}_{+}+{}_{+}=-\frac{}{h}& \left(15\right)\end{array}$

(98) Where non-zero entries of the diagonal matrices and may be computed as:

(99) $\begin{array}{cc}{}_{i,i}=\frac{1}{\left(1+h^{-1}{C}_{i,i}\right)B_{i,i}^{}}& \left(16\right)\\ {}_{i,i}=\frac{h^{-2}{C}_{i,i}}{\left(1+h^{-1}{C}_{i,i}\right)B_{i,i}^{}}& \left(17\right)\end{array}$

(100) This may be related to how constrained multibody simulators may allow combination of compliant constraints and Baumgarte stabilization to be interpreted as an implicit stiff spring and damper.

(101) Turning now to FIG. 3, a simulation 302 of a heavy ball (i.e., a 100 Kg load) attached to a flexible cable being dropped and beginning swinging with a pendulum motion is illustrated. A first graph 304 shows kinetic energy when simulated using the method of Tournier et al. according to a first plot 306 and using the method of the present technology according to a second plot 308.

(102) Turning now to FIG. 4, a second graph 402 and a third graph 404 illustrate instantaneous work at each frame using an original geometric stiffness matrix and a diagonal approximation of the original geometric stiffness matrix to simulate the heavy ball of FIG. 3 using hinge joints (the second graph 402) and a flexible cable joint (the third graph 404). FIG. 4 demonstrates that the approximation set forth above in connection with the equation (12) bounds the instantaneous work. For rotational joints, such as a hinge, the diagonal approximation tends to match the instantaneous work relatively closely, whereas for examples modeled using a flexible cable joint, the instantaneous work may be somehow overestimated but nonetheless achieves the goal of upper bounding the work done by {tilde over (K)}.

(103) Turning now to FIG. 5, a fourth graph illustrating performance gains when using low rank updates to compute {tilde over (M)}.sup.1. As it may be appreciated, for angular joints, the gain may appear to scale linearly with the number of joints, in accordance with table 2 below illustrating an average time to compute {tilde over (M)}.sup.1 for a cable using various numbers of ball-and-socket joints.

(104) TABLE-US-00002 TABLE 2 Number of joints Method 10 50 100 Full rank 0.099 ms 6.795 ms 49.668 ms Low rank 0.010 ms 0.517 ms 2.633 ms Gain 9.81 13.13 18.87

(105) Turning now to FIG. 6, a scene 600 represents a character with super strength 602 preventing a tank 604 from driving away. The scene 600 involves many different joints, including hinges, ball-and-sockets, prismatic and flexible cable joints. FIG. 6 also comprises a fifth graph 606 illustrating performance gains using a low rank inverse for the scene 600. As it may be appreciated, the geometric stiffness matrix for the cable joints is not sparse and, as a result, performance gains may be diminished as a number of joints is increased.

(106) Turning now to FIG. 7, a sixth graph 702 illustrates performance gain and computation time for cables simulated using various numbers of flexible cable joints. A speedup from low rank updates is mitigated due to a non-sparse structure of the geometric stiffness matrix.

(107) Turning now to FIG. 8, a cable attached to a heavy load is represented in multiple states 802. The states 802 illustrate vibrations in the cable when integrated with a small time step and no damping h= 1/100,000 s (left) and large time step h=0.01 s with adaptive damping scheme (middle) and without (right). As may be seen from FIG. 8, dynamical behaviour may be better preserved with the adaptive damping scheme. FIG. 8 also illustrates a seventh graph 804 which illustrates damping coefficients of a selected link for the cable simulation of FIG. 3. For less massive objects (10 Kg) stabilization damping is zero for most of the simulation. As the mass increase, more damping is required as constraint forces increase.

(108) Turning now to FIG. 9, scenes 902, 904 and 906 illustrate an overlay of selected frames from a crane simulation in accordance with the present technology. A 100,000 Kg boom is elevated 20 degrees and then dropped. The simulation remains stable for 0.1 s (902), 0.01 s (904) and 0.001 s (906) time steps.

(109) FIGS. 10 and 11 illustrate scenes 1002, 1004, 1006, 1008 and 1010 generated by a simulation in accordance with the present technology. An heavy vehicle attached to a nearly inextensible cable is dropped from an end of a crane. The simulation involves very large mass ratios and heterogeneous collection of joints. The simulation remains stable at a time step of 1/60 s. As it may be appreciated, cable dynamics may be better preserved by the adaptive damping scheme of the present technology.

(110) Results and Discussions

(111) This section provides examples of evaluation of some aspects of the present technology. Some examples aim at examining some of the computational gains from using a method of low rank matrix inversion in accordance with the present technology whereas other examples illustrate potential benefits of some embodiments of the adaptive damping approach of the present technology. Results were obtained using an Intel Core i73.3 GHz CPI from Intel and a single thread. Unless otherwise stated, time steps used for experiments is h=0.01 s. The physics engine Vortex from CM Labs has been used for the tests. A linear in accordance with the equation (3) is used for all simulations and solved using a standard Cholesky decomposition. The same linear system may also be solved by LU decomposition.

(112) Examples illustrated in FIGS. 3 and 8 use a constraint space damping method whereas as other examples are stabilized by inertial damping.

(113) Performance of Low Rank Updates

(114) Three examples are discussed in the paragraphs below to evaluate computational speedup that may be obtained by using low rank updates in accordance with the present technology.

(115) Articulated Chains

(116) A computational time of low rank updates and full rank matrix inversion is compared in FIG. 5. Articulated chains involving various numbers of joints and joints types are simulated and a computational time required to compute {tilde over (M)}.sup.1 is measured. As indicated by table 2, performance gain grows as a number of joints in the system increases. Likewise, a gain is consistent across hinge, universal and ball-and-socket joints.

(117) Strong Robot

(118) An example involving a heterogeneous collection of joints is shown in FIG. 6. The robot character is modeled using compliant constraints, including 5 ball-and-sockets, 11 universal and 10 hinge joints. Each foot of the robot is anchored to the ground using a hinge joint and each fingertip is attached to an end of the cable using a ball-and-socket joint. Masses range from 0.1 Kg for the finger bodies to 15,000 Kg for the vehicle attached to an opposite end of the cable. The vehicle is driven by a motorized prismatic joint to move it away from the robot. The DOFs in the robot are controlled using implicit PD servos with a set point matching the original pose of the character. This scenario resembles a character animation sequence that might occur in a video game with many different types of joints and large mass ratios. Nevertheless, the simulation remains stable and the motion of the character is smooth.

(119) Flexible Cable Joint

(120) A special joint may be used to model the cable for the scene depicted in FIG. 6. The joint uses several dot-2 constraints, for which the geometric stiffness matrix may be relatively dense. The graph 606 of FIG. 6 shows that performance gains diminish as the number of flexible joints is increased and eventually results in longer computation times. Performance results from these simulations are shown in FIG. 7.

(121) Validation

(122) Evaluation of an adaptive damping method in accordance with the present technology may be conducted by monitoring energy of the system, vibrations in chains and explore the range of damping coefficients for systems with various mass ratios.

(123) Kinetic Energy

(124) Simulation of a pendulum swinging under gravity exchanges potential energy and kinetic energy. If friction and other dissipative elements in the simulation are nominal, the total mechanical energy should be nearly constant. FIG. 3 shows the kinetic energy per time step of a cable with 100 Kg load and pendulum motion. The flexible cable joints have an axial stiffness of 10.sup.12 and bending and torsional stiffness of 10. Using known methods, such as Tournier et al., dissipation caused by geometric stiffness causes a peak kinetic energy to decrease by nearly 40% after 40 seconds of simulation time. When the stability criterion with constraint damping is used, the peak kinetic energy drops by only 2% after the same amount of time, indicating that the mechanical energy in the system is better preserved.

(125) Cable Vibrations

(126) In FIG. 8, vibrations of a cable are visualized when an attached heavy load (in this case 5000 Kg) is dropped straight down and reaches its apex. A simulated behavior with very small time steps and no geometric stiffness or damping is provided for comparison. The vibrations seen when the adaptive damping method in accordance with the present technology is used better represent the observed oscillations in a high fidelity simulation and without the dynamical behavior appears to be over damped.

(127) Damping Coefficients

(128) The graph 804 of FIG. 8 shows the damping coefficients of a selected cable link for the example shown in FIG. 3. The mass is increased from 10 Kg to 100 Kg to 1000 Kg. This may demonstrate that the adaptive damping in accordance with the present technology adapts to constraint forces at a current simulation frame. For smaller masses, damping is only required when load reaches its apex (around frame 160) and constraint forces are largest. As the load increases, damping may be required more often to stabilize the system

(129) Other results discussed below demonstrate that the present technology may also be effective by producing stable and interactive simulations with large mass ratios and/or large time steps.

(130) Crane With Heavy Load

(131) FIGS. 10 and 11 shows a 15,000 Kg vehicle attached to a cable being dropped from a crane. The cable is nearly inextensible, with an axial stiffness of 10.sup.10. The mass ratio in the system is similarly large, ranging from 100,000 Kg for the boom to 0.2 Kg for individual cable segments. Despite these challenging parameters, the simulation is stable when integrated using a time step of 1/60 s. The adaptive damping method according to the present technology may help to preserve rich cable dynamics and may be applied to angular velocities of simulation bodies with =1.0. The articulated linkage of the crane is modeled using hinge joints, whereas the cable uses relaxed rigid joints.

(132) Varying Time Steps

(133) A damping coefficient computed in accordance with the present technology may be dependent on the constraint forces and masses of bodies in the simulation. However, the equation (11) may suggest that the damping coefficient may also be dependent on the time step. FIG. 9 illustrates motion of the crane boon when it is dropped from a 20 degrees angle of elevation. Despite large mass ratios and stiff constraints, the simulation remains stable at time steps of 0.1 s, 0.01 s and 0.001 s. Also, the dynamical behavior in each case may be similar, although for the largest time step value the simulation may not appear to be more damped and may exhibit less high frequency dynamics.

(134) Constraint Space Vs Inertial Damping

(135) For cable simulations, artifacts may be observed when using inertial damping formulation. Specifically, there may be dissipation about the torsional axis.

(136) The constraint space damping method in accordance with the present technology may be well suited for 6D rigid and flexible cable joints. In some embodiments, constraint space may be augmented by additional constraint rows with negligible compliance values. This may be done automatically at each time step or manually when the simulation model is being designed.

(137) In some embodiments, articulated simulations may be stabilized by only damping rotational DOFs. For many applications, such as physics-based character animation only hinges, universal and/or ball-and-socket joints may be used to model characters. We may however observe that for mechanisms using the dot-2 constraint, such as the prismatic and the flexible cable joint, applying damping to the translational degrees of freedom may result in a behavior where bodies appear to move through a viscous fluid. For such cases, we may simply discard the translational damping coefficient and we may set B.sub.i,i=0 for these DOFs. This may help to prevent viscosity artifacts, especially for cable simulation, but may have similar stabilizing properties.

(138) As it may be appreciated, some aspects of the present technology may provide at least some benefits. Using the geometric stiffness in its diagonal form, a stability criterion may be used to adaptively introduce a minimal amount of damping a three fold may benefit: (i) the system matrix may remain positive definite, permitting more efficient numerical solvers; (ii) the effort of manually tuning artificial damping coefficients may be reduced and (iii) stable simulation systems with unprecedented stiffness and/or mass rations may be possible at real-time frame rates.

(139) Constraint Library

(140) This section contains derivations of the geometric stiffness for a basic set of constraints that may be commonly used in articulated rigid body simulation. Each derivation considers the case of a two body system where the geometric stiffness may be represented by a 1212 matrix. For succinctness, the 33 block structure of theses matrices is exploited.

(141) Basic Partial Derivatives

(142) The instantaneous angular velocity of a body and the time derivative of its rotation matrix {dot over (R)} are related by
{circumflex over ()}={dot over (R)}R.sup.1.(18)

(143) We can reinterpret this lemma for small angular displacements such that
{circumflex over ()}=dRR.sup.1,
and post-multiplying by R gives
{circumflex over ()}R=dR.(19)

(144) Furthermore, for a vector n.sup.3 attached to a rotating body such that n=Rn.sup.1, the spatial derivative of the vector is given by

(145) $\begin{array}{cc}\frac{n}{}=\hat{n}.& \left(20\right)\end{array}$

(146) Equations (19) and (20) may be the building blocks used to derive the geometric stiffness equations for a basic constraint library. Additional mathematical details of rigid body kinematics may be found in Murray et al. (MURRAY R. M., LI Z., SASTRY S. S., SASTRY S. S.: A mathematical introduction to robotic manipulation. CRC press, 1994.), the entirety of which is hereby incorporated by reference in jurisdictions allowing such incorporation.

(147) Ball- and Socket Joint

(148) A ball-and-socket, or spherical, joint connecting bodies i and j constrains two points on the bodies to have the same position while allowing relative angular motion. The constraint equation can be written as
.sub.bs=p.sub.i+s.sub.ip.sub.js.sub.j=0.

(149) Rearranging Eq. 18 gives {dot over (R)}={circumflex over ()}R, and so She time derivative of the constraint equation can be written as
.sub.bs={dot over (p)}.sub.i.sub.i.sub.i{dot over (p)}.sub.j+.sub.j.sub.j,
giving the constraint Jacobian matrix
J.sub.bs=(I.sub.iI .sub.j).

(150) Its other words,

(151) ${\stackrel{.}{}}_{\mathrm{bs}}=J_{\mathrm{bs}}\left(\begin{array}{c}{\stackrel{.}{p}}_i\\ {}_i\\ {\stackrel{.}{p}}_j\\ {}_j\end{array}\right).$

(152) The constraint forces acting on fee two bodies are

(153) $\begin{array}{cc}{J}_{\mathrm{bs}}^T{}_{\mathrm{bs}}=\left(\begin{array}{c}{}_{\mathrm{bs}}\\ {\hat{s}}_i{}_{\mathrm{bs}}\\ -{}_{\mathrm{bs}}\\ -{\hat{s}}_i{}_{\mathrm{bs}}\end{array}\right),& \left(21\right)\end{array}$
where .sub.bs .sup.3 are the Lagrange multipliers for the ball-and-socket constraint. Differentiating Eq. 21 with respect to the spatial configuration of the bodies and applying the chain rule gives

(154) $\frac{J_{\mathrm{bs}}^T{}_{\mathrm{bs}}}{x}=J_{\mathrm{bs}}^T\frac{{}_{\mathrm{bs}}}{x}+\frac{J_{\mathrm{bs}}^T}{x}{}_{\mathrm{bs}},$
where the first term is zero and the second term is the geometric stiffness matrix of the joint, or

(155) 0${\tilde{K}}_{\mathrm{bs}}=\frac{J_{\mathrm{bs}}^T}{x}{}_{\mathrm{bs}}.$
Expanding this into 33 blocks gives

(156) $\begin{array}{cc}{\tilde{K}}_{\mathrm{bs}}=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& {\hat{}}_{\mathrm{bs}}{\hat{s}}_i& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& -{\hat{}}_{\mathrm{bs}}{\hat{s}}_j\end{array}\right).& \left(22\right)\end{array}$

(157) Noting that s.sub.i=R.sub.is.sub.i and s.sub.j=R.sub.js.sub.j, equation (22) matches the matrix provided in Tournier et al.

(158) Axial Lock

(159) An axial lock constrains the angular motion about a single axis. The constraint Jacobian matrix is the 12 component row vector
J.sub.n(0n.sup.T.sub.i0n.sup.T.sub.i)
where the locked axis of rotation n.sub.i is fixed to body i. The geometric stiffness for this constraint is

(160) $\begin{array}{cc}{\tilde{K}}_{\mathrm{ax},n}=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& -{}_{\mathrm{ax}}{\hat{n}}_i& 0& 0\\ 0& 0& 0& 0\\ 0& {}_{\mathrm{ax}}{\hat{n}}_i& 0& 0\end{array}\right).& \left(23\right)\end{array}$

(161) Universal Joint

(162) The universal, or Hooke, joint has similar behavior to the spherical joint, but removes a rotational degree of freedom by keeping an axis fixed to body i perpendicular fixed to body j. The constraint equations for this joint can be written as
.sub.bs=p.sub.is.sub.ip.sub.js.sub.j=0
.sub.dl,u=n.sup.T.sub.iu.sub.j=0
where n.sub.i is a unit vector in coordinate frame of body i which should remain orthogonal to unit vector u.sub.j in coordinate frame of body j. Since the derivation of .sub.bs is the same as previously discussed, we instead focus on the dot-1 constraint .sub.d1. The Jacobian matrix of the dot-1 constraint is
J.sub.d1,u=(0({circumflex over (n)}.sub.iu.sub.j).sup.T0({circumflex over (n)}.sub.iu.sub.j).sup.T).

(163) Note that for vectors a, b.sup.3 that b={circumflex over (b)}a and (b).sup.T=b.sup.T. The geometric stiffness for the dot-1 constraint is

(164) $\begin{array}{cc}{\tilde{K}}_{d1,u}=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& & 0& -{}^T\\ 0& 0& 0& 0\\ 0& -& 0& {}^T\end{array}\right)& \left(25\right)\end{array}$
where =.sub.d1.sub.j.sub.i and the scalar value .sub.d1 is the Lagrange multiplier of the dot-1 constraint at Else previous time step. The geometric stiffness for the universal joint may be assembled as
K.sub.un=K.sub.bs+K.sub.d1,u.

(165) Hinge Joint

(166) The hinge, or revolute, joint allows a relative rotation of two bodies about a single axis. The joint behaves much like a universal, but with an additional dot-1 constraint equation
.sub.bs=p.sub.i+s.sub.ip.sub.js.sub.j=0
.sub.d1,u=n.sub.i.sup.Tu.sub.j=0
.sub.d1,v=n.sub.i.sup.Tv.sub.j=0
where v.sub.j is a unit vector in coordinate frame of body j which is orthogonal to u.sub.j and is constrained to be perpendicular to n.sub.i. The geometric stiffness {tilde over (K)}.sub.d1,v for the dot-1 constraint between vector v.sub.j and n.sub.i is found by applying Eq. 25 and substituting v.sub.j or u.sub.j. The hinge geometric stiffness is then assembled as
{tilde over (K)}.sub.ni={tilde over (K)}.sub.bs+{tilde over (K)}.sub.d1,u+{tilde over (K)}.sub.d1,v.

(167) Prismatic Joint

(168) A prismatic joint allows a relative translation between two bodies along a single axis. The constraint equations for this joint can be written as
.sub.d1,u=n.sub.i.sup.Tu.sub.j=0
.sub.d1,v=n.sub.i.sup.Tv.sub.j=0
.sub.d1,n=u.sub.i.sup.Tn.sub.j=0
.sub.d2,u=d.sub.ij.sup.Tu.sub.i=0
.sub.d2,v=d.sub.ij.sup.Tv.sub.i=0
where d.sub.ij=p.sub.i+s.sub.ip.sub.js.sub.j. This joint consists of three dot-1 constraints (top rows), plus equations .sub.d2,u and .sub.d2,v, which are dot-2 constraints, A single Lagrange multiplier is used for each dot-2 constraint, which together eliminate relative translation of the bodies along any axes perpendicular to n.sub.i.

(169) Using .sub.d2,u as an example, which constrains relative motion in the direction u.sub.i, the Jacobian matrix is
J.sub.d2,u=(u.sub.i.sup.Tu.sub.i.sup.T({circumflex over (d)}.sub.ij.sub.i)u.sub.i.sup.Tu.sub.i.sup.T.sub.j)(26)
with the constraint Jacobian for .sub.d2,v having a similar form. The geometric stiffness matrix fat J.sub.d2,u is

(170) $\begin{array}{cc}{\tilde{K}}_{d2,u}={}_{d2,u}\left(\begin{array}{cccc}0& {\hat{u}}_i^T& 0& 0\\ {\hat{u}}_i& \left({\hat{d}}_{i,j}-{\hat{s}}_i\right){\hat{u}}_i& {\hat{u}}_i^T& {\left({\hat{s}}_j{\hat{u}}_i\right)}^T\\ 0& {\hat{u}}_i& 0& 0\\ 0& {\hat{s}}_j{\hat{u}}_i& 0& -{\left({\hat{s}}_j{\hat{u}}_i\right)}^T\end{array}\right).& \left(27\right)\end{array}$

(171) The matrix {tilde over (K)}.sub.d2,v may be found by applying Eq. 27 and substituting v.sub.i for u.sub.i. The geometric matrix for the prismatic joint is competed as
{tilde over (K)}.sub.pr={tilde over (K)}.sub.d1,u+{tilde over (K)}.sub.d1,v+{tilde over (K)}.sub.d1,n+{tilde over (K)}.sub.d2,u+{tilde over (K)}.sub.d2,v.

(172) Rigid Joint

(173) The prismatic joint is augmented with a sixth constrains equation that restricts the translational motion of the bodies, effectively constraining the rotational and translational motion of the two bodies to be file same. This is done using an additional dot-2 constraint.
.sub.d2,n=d.sub.ij.sup.Tn.sub.i=0.

(174) The geometric stiffness matrix of the rigid joint is then computed as
{tilde over (K)}.sub.6D={tilde over (K)}.sub.pr+{tilde over (K)}.sub.d2,u.

(175) Flexible Cable Joint

(176) By choosing small or moderate compliances for the rigid joint, the constraint rows are relaxed and the joint becomes flexible. Compliance and damping values are then selected for each of the three rotational and three translational constrained degrees of freedom to match the behavior of an extensible cable. For example, Eq. 16 and Eq. 17 are tuned in this way for the cable results shown in Section 5.

(177) Low Rank Decompositions

(178) This section provides low rank decompositions of basic constraint library set forth above. All factorizations consider the case of a two body system with 1212 geometric stiffness matrix, although decompositions for larger systems may be possible.

(179) Ball- and Socket

(180) The low rank factorization of {tilde over (K)}.sub.bs follows from the analysis in Section 3.1. The vectors {right arrow over (r)}, {right arrow over (s)}, {right arrow over (t)} and the Lagrange multipliers .sub.r, .sub.s, .sub.t are used with subscripts to indicate the body index. Using this notation, the decomposition is as follows:
s.sub.i=(0 {right arrow over (s)}.sub.i.sup.T 0 0).sup.T
r.sub.i=(0 {right arrow over (r)}.sub.i.sup.T 0 0).sup.T
t.sub.i=(0 {right arrow over (t)}.sub.i.sup.T 0 0).sup.T
s.sub.j=(0 0 0 {right arrow over (s)}.sub.j.sup.T).sup.T
r.sub.j=(0 0 0 {right arrow over (r)}.sub.j.sup.T).sup.T
t.sub.j=(0 0 0 {right arrow over (t)}.sub.j.sup.T).sup.T
{tilde over (K)}.sub.bs=s.sub.i(.sub.r.sub.is.sub.ir.sub.i.sup.T.sub.s.sub.i{right arrow over (r)}.sub.i{right arrow over (r)}.sub.i.sup.T+.sub.r.sub.is.sub.it.sub.i.sup.T.sub.s.sub.it.sub.it.sub.i.sup.T)s.sub.j(.sub.r.sub.js.sub.jr.sub.j.sup.T.sub.s.sub.j{right arrow over (r)}.sub.j{right arrow over (r)}.sub.j.sup.T+.sub.r.sub.j s.sub.jt.sub.j.sup.T.sub.s.sub.jt.sub.jt.sub.j.sup.T).

(181) Axial Constraint

(182) The axial constraint decomposition has the following form:
.sub.i=(0 u.sub.i.sup.T 0 0).sup.T
v.sub.i=(0 v.sub.i.sup.T 0 0).sup.T
.sub.i=(0 0 0 u.sub.i.sup.T).sup.T
v.sub.i=(0 0 0 v.sub.i.sup.T).sup.T
{tilde over (K)}.sub.ax=(.sub.iv.sub.i.sup.T+v.sub.i.sub.i.sup.T+.sub.iv.sub.i.sup.Tv.sub.i.sub.i.sup.T).
Dot-1 Constraint

(183) The decomposition for the dot-1 constraint follows the Jacobian definition from Eq. 24, using vectors n.sub.i and u.sub.j, but the decomposition is valid for other vector pairings. The low rank decomposition of the dot-1 geometric stiffness matrix is:
n=(0 0 0n.sub.i.sup.T).sup.T
=(0 u.sub.j.sup.T 0 0).sup.T
n=(0 n.sub.i.sup.T 0 n.sub.i.sup.T).sup.T
=(0 u.sub.j.sup.T 0 u.sub.j.sup.T).sup.T
{tilde over (K)}.sub.d1,u=(n.sup.T+n.sup.T).

(184) Dot-2 Constraint

(185) Finally, for the dot-2 constraint we begin by forming an orthonormal basis using the vector

(186) ${\stackrel{.fwdarw.}{d}}_{\mathrm{ij}}=\frac{d_{\mathrm{ij}}}{||d_{\mathrm{ij}}||},$
such that
{right arrow over (d)}.sub.ij{right arrow over (e)}.sub.ij,{right arrow over (e)}.sub.ij{right arrow over (f)}.sub.ij,{right arrow over (d)}.sub.ij{right arrow over (f)}.sub.ij, {right arrow over (d)}.sub.ij={right arrow over (e)}.sub.ij={right arrow over (f)}.sub.ij=1.

(187) This basis, along with the orthonormal bases {right arrow over (r)}.sub.i, {right arrow over (s)}.sub.i, {right arrow over (t)}.sub.i and {right arrow over (r)}.sub.j, {right arrow over (s)}.sub.j, {right arrow over (t)}.sub.j, is projected onto the constraint coordinate vectors v.sub.i, and n.sub.i giving the coefficients

(188) TABLE-US-00003 .sub.i, n = {right arrow over (r)}.sub.i.sup.T n.sub.i .sub.i, v = {right arrow over (r)}.sub.i.sup.T v.sub.i .sub.i, n = {right arrow over (t)}.sub.i.sup.T n.sub.i .sub.i, v = {right arrow over (t)}.sub.i.sup.T v.sub.i .sub.j, n = {right arrow over (r)}.sub.j.sup.T n.sub.i .sub.j, v = {right arrow over (r)}.sub.j.sup.T v.sub.i .sub.j, n = {right arrow over (t)}.sub.j.sup.T n.sub.i .sub.j, v = {right arrow over (t)}.sub.j .sup.T v.sub.i .sub.n = {right arrow over (f)}.sub.ij.sup.Tn.sub.i .sub.v = {right arrow over (f)}.sub.ij.sup.Tv.sub.i .sub.n = {right arrow over (e)}.sub.ij.sup.Tn.sub.i .sub.v = {right arrow over (e)}.sub.ij.sup.Tv.sub.i.

(189) We also define the following vectors:

(190) TABLE-US-00004 n.sub.2 = (0 n.sub.i.sup.T 0 0).sup.T v.sub.2 = (0 v.sub.i.sup.T 0 0).sup.T n.sub.31 = (n.sub.i.sup.T 0 n.sub.i.sup.T 0).sup.T v.sub.13 = (v.sub.i.sup.T 0 v.sub.i.sup.T 0).sup.T n.sub.42 = (0 n.sub.i.sup.T 0 n.sub.i.sup.T).sup.T v.sub.24 = (0 v.sub.i.sup.T 0 v.sub.i.sup.T).sup.T r.sub.i = (0 {right arrow over (r)}.sub.i.sup.T 0 0).sup.T t.sub.i = (0 {right arrow over (t)}.sub.i.sup.T 0 0 ).sup.T r.sub.j = (0 0 0 {right arrow over (r)}.sub.j.sup.T).sup.T t.sub.j = (0 0 0 {right arrow over (t)}.sub.j.sup.T).sup.T f.sub.ij = (0 {right arrow over (f)}.sub.ij.sup.T 0 0).sup.T .sub.ij = (0 {right arrow over (e)}.sub.ij.sup.T 0 0).sup.T c.sub.1 = .sub.j, nt.sub.j + .sub.j, n r.sub.j d.sub.1 = .sub.j, vt.sub.j ++ .sub.j, vr.sub.j c.sub.2 = .sub.j, n v.sub.2 .sub.j, vn.sub.2 d.sub.2 = .sub.j, n v.sub.2 .sub.j, vn.sub.2 c.sub.4 = .sub.n v.sub.2 .sub.vn.sub.2 d.sub.4 = .sub.n v.sub.2 .sub.vn.sub.2 c.sub.5 = .sub.i, n v.sub.2 .sub.i, vn.sub.2 d.sub.5 = .sub.i, n v.sub.2 .sub.i, vn.sub.2.

(191) The low rank decomposition {tilde over (K)}.sub.d2,u becomes
{tilde over (K)}.sub.d2,u=(v.sub.13n.sub.2.sup.T+n.sub.31v.sub.2.sup.T+n.sub.2v.sub.13.sup.T+v.sub.2n.sub.31.sup.T)+s.sub.j(n.sub.42d.sub.1.sup.T+v.sub.24c.sub.1.sup.T+r.sub.jd.sub.2.sup.Tt.sub.jc.sub.2.sup.T)+d.sub.ij(f.sub.ijd.sub.4.sup.T.sub.ijc.sub.4.sup.T)s.sub.i(r.sub.id.sub.5.sup.Tt.sub.ic.sub.5.sup.T).

(192) Turning now to FIG. 12, a flowchart illustrating a computer-implemented method 1200 of executing a simulation of a constrained multi-body system is illustrated. In some embodiments, the computer-implemented method 1200 may be (completely or partially) implemented on a computing environment similar to the computing environment 100.

(193) The method 1200 comprises using a physics engine simulating the constrained multi-body system. The constrained multi-body system comprises articulated constraints, the articulated constraints are associated with a geometric stiffness matrix. The geometric stiffness matrix defines a geometric stiffness. At a step 1202, the method 1200 executes generating a diagonal approximation of the geometric stiffness matrix. Then, at a step 1204, the diagonal approximation is used as part of a stability analysis in which damping is automatically adjusted so that the damping stabilizes the simulation of the constrained multi-body system.

(194) In some embodiments, the damping is associated with a damping coefficient. The damping coefficient may be automatically adjusted so as to be large enough to stabilize transverse oscillations of the simulation of the constrained multi-body system. In some embodiments, the damping is an adaptive damping, the adaptive damping being associated with damping values calculated at each time step. The damping values may be calculated based on the following equation:

(195) $b=\frac{h^{2}k-4m}{2h}.$

(196) Wherein k is a stiffness, b is a damping and m is a mass, h is a given time step and a is a positive scalar value that allows control over how close dynamical system may be to the boundary before damping may be required.

(197) In some embodiments, the simulation of the constrained multi-body system comprises a dynamical system wherein the geometric stiffness matrix is reinterpreted as a spring and the generated forces are used in combination with a stability criterion for stable simulation. The spring forces may define the damping.

(198) In some embodiments, the stability criterion may comprise a threshold based on a behavior of a simple 1D undamped oscillator. The geometric stiffness matrix may serve as the stability criterion as to whether the dynamical system becomes unstable. In some embodiments, the geometric stiffness matrix is a {tilde over (K)} matrix, wherein the {tilde over (K)} matrix is defined as follows:

(199) $\tilde{K}=\frac{J^T}{x}.$

(200) Wherein J is a matrix of constraint Jacobians, are constraint forces and x are positions and orientations of bodies in the simulation of the constrained multi-body system.

(201) In some embodiments, the diagonal approximation of the geometric stiffness matrix may preserve mechanical work properties associated with the original geometric stiffness matrix. The geometric stiffness may define how forces in the constrained multi-body system change due to instantaneous motion.

(202) In some embodiments, the instantaneous motion comprises at least one of rotation of links associated with the constrained multi-body system, translation of links associated with the constrained multi-body system and relative motion of links associated with the constrained multi-body system.

(203) In some embodiments, a topology of a mechanical structure of the constrained multi-body system may remain unmodified during the simulating of the constrained multi-body system. In some embodiments, the parameters of the damping define an interpretation of the geometric stiffness used in the simulation of the constrained multi-body system. In some embodiments, at least one of the parameters define physical behaviour of a transient physical spring that is able to generate forces in directions transverse to transverse oscillations directions. In some embodiments, the simulating of the constrained multi-body system may be executed in real-time.

(204) In some embodiments, stabilizing the simulation of the constrained multi-body system may comprise stabilizing the transverse oscillations of the simulation of the constrained multi-body system.

(205) In some embodiments, the geometric stiffness matrix may be generated by using a low rank update algorithm, the using of the low rank update algorithm comprising constructing a lead matrix of dynamical equations of motion, the lead matrix comprising the geometric stiffness matrix. The geometric stiffness matrix may model the articulated constraints. In some embodiments, data encoding parameters of the transient spring models a tensor encoding variations in constraint force directions.

(206) As it may be appreciated, a constrained multi-body system may comprise rigid bodies such as the ones illustrated in FIG. 2-11. In some embodiments, a rigid body may be defined as an idealization of a solid body in which deformation is neglected. In other words, a distance between any two given points of a rigid body may remain constant in time regardless of external forces exerted on it. In a constrained multi-body system simulation, constraints may specify how the rigid bodies are allowed to move relative to one another.

(207) As it may be appreciated, a physics engine may compute motions of rigid bodies based on dynamic equations, typically subject to gravity and possibly connected using kinematic constraints also referred to as joints, such as, but not limited to, a hinge or slider constraint. Examples of physics engine may include, but is not limited to, Bullet, Open Dynamics Engine (ODE), Havok, Simmechanics from Matlab and Vortex from CM Labs.

(208) As it may be appreciated, an articulated constraint may consist of a constraint function defining allowed relative motion between rigid bodies. A derivative of the constraint may be a Jacobian matrix. The derivative of the Jacobian w.r.t. position may be the geometric stiffness matrix in accordance with the present technology.

(209) In some embodiments, forces described by the geometric stiffness may be considered as forces generated by a set of springs in relevant dimensions. Dynamical equations of motion may be rewritten in such a way as to match equations of a spring, thus providing stiffness and damping coefficient of such a spring.

(210) In some embodiments, stability may depend on multiple factors such as mass, stiffness, damping, the time step, and the integration scheme. The stability analysis may assume that a system is linear and may use standard methods to determine a boundary of stability. The stability analysis may be based on the fact that if the magnitude of the Eigen values of the phase space matrix are <1 (for an example of such a matrix, see P from equation (8)). We may do this for a simple harmonic oscillator, and for such a simple system the stability may be guaranteed if this condition holds. The present technology may also be applied to more complex systems using an approximation of the geometric stiffness thereby enlarging the region stability.

(211) In some embodiments, the stiffness matrix may be diagonalized to increase performance of the updated mass matrix computation and facilitate convergence analysis. In some embodiment, the implicit damping matrix may be then computed using equation (13).

(212) In some embodiments, damping, adaptive damping and adaptive damping scheme may refer to a method for varying an amount of damping in a damper according to some criterion. In some embodiments, an amount of damping may be different for each degree of freedom and may depend on an amount of tension in system in that degree of freedom, and be only as high as necessary to maintain stability. In some embodiments, the method for varying an amount of damping may mitigate non-constitutive energy dissipation.

(213) While the above-described implementations have been described and shown with reference to particular steps performed in a particular order, it will be understood that these steps may be combined, sub-divided, or re-ordered without departing from the teachings of the present technology. At least some of the steps may be executed in parallel or in series. Accordingly, the order and grouping of the steps is not a limitation of the present technology.

(214) It should be expressly understood that not all technical effects mentioned herein need to be enjoyed in each and every embodiment of the present technology.

(215) Modifications and improvements to the above-described implementations of the present technology may become apparent to those skilled in the art. The foregoing description is intended to be exemplary rather than limiting. The scope of the present technology is therefore intended to be limited solely by the scope of the appended claims.