{"id":150,"date":"2015-10-02T16:05:28","date_gmt":"2015-10-02T16:05:28","guid":{"rendered":"http:\/\/project.inria.fr\/softrobot\/?page_id=150"},"modified":"2017-04-03T15:07:41","modified_gmt":"2017-04-03T15:07:41","slug":"real-time-inverse-simulation-method","status":"publish","type":"page","link":"https:\/\/project.inria.fr\/softrobot\/documentation\/real-time-inverse-simulation-method\/","title":{"rendered":"Real Time Inverse Simulation Method"},"content":{"rendered":"

Introduction<\/strong><\/h1>\n

This section details the formulation that is used for real-time inverse method on FEM. It describes how it can\u00a0be used to estimate the external loads, pressures, displacements… that lead to a given\u00a0deformation. The static FEM accounts for non-linear geometrical deformations\u00a0and integrates over the structure a Hookean constitutive law. During\u00a0each step i of the simulation, a linearization of the internal forces is computed:<\/p>\n

f(xi<\/sub>) \u2248 f(xi\u22121<\/sub>) + K(xi\u22121<\/sub>)dx \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 \u00a0(1)<\/p>\n

where f provides the volumetric internal stiffness forces at a given position x\u00a0of the nodes, K(x) is the tangent stiffness matrix that depends on the actual position of the nodes and dx is the difference between consecutive positions in time dx = xi<\/sub>\u00a0\u2212 xi\u22121<\/sub>. The lines and columns that correspond to fixed nodes are removed from the system to improve the condition number of the matrix K. Static equilibrium (the sum of external and internal forces equal to zero) is sought at each step:<\/p>\n

\u2212K(xi\u22121<\/sub>)dx = p + f(xi\u22121<\/sub>) + JT<\/sup>\u00a0\u03bb \u00a0 \u00a0 \u00a0 \u00a0(2)<\/p>\n

where p accounts for the external forces (e.g. gravity) that are known and JT<\/sup>\u00a0\u03bb gathers the\u00a0contributions of the Lagrange multipliers. Two types of multipliers are defined:<\/p>\n

actuators multipliers \u03bba<\/sub><\/strong>:\u00a0we use these constraints to describe the external efforts applied on the boundaries (such as pressure for instance). The\u00a0location of the boundary conditions are supposed to be known, the directions of\u00a0the effort JT\u00a0<\/sup>can be updated at each step, and \u03bba<\/sub>\u00a0is the unknown intensity of\u00a0the efforts. We can set (and update at each step i) an interval of\u00a0potential values as min \u2264 \u03bba<\/sub>\u00a0\u2264 max.<\/p>\n

effectors multipliers \u03bbe<\/sub><\/strong>: as we don’t want to act on the effectors but only on the actuators, we assign \u03bbe<\/sub>\u00a0= 0.\u00a0Even if they are null, these multipliers are useful to\u00a0build the optimization problem.<\/p>\n

Indeed, the next step consists in the projection of the FEM model equations\u00a0into the constraint space: the size of matrix K is often very large so an optimization\u00a0in the motion space would be very computationally expensive. Instead,\u00a0using the Schur complement of the constraint problem, we do a projection that\u00a0drastically reduces the size of the search space.\u00a0Thus three steps are performed.<\/p>\n

Step I:<\/strong>\u00a0a free configuration xfree<\/sup>\u00a0of the deformable model is found by solving\u00a0equation (2) with \u03bb = 0. For the constraints defined on effectors, we\u00a0compute a violation noted \u03b4e<\/sub>free<\/sup>. For instance, in a position effector,\u00a0\u03b4e<\/sub>free<\/sup>\u00a0provides a vector between the effector position during the free motion and the desired position.<\/p>\n

Step II<\/strong>: This step is central in the method. It consists in projecting the mechanics into the constraint space. As the constraints are the inputs (effectors) and outputs (actuators) of the inverse problem, we obtain the smallest possible projection space for the inverse problem:<\/p>\n

\u03b4e<\/sub>\u00a0= Je<\/sub>K\u22121<\/sup>Ja<\/sub>T<\/sup>\u03bba<\/sub>\u00a0 + \u03b4e<\/sub>free<\/sup>\u00a0 \u00a0 \u00a0 \u00a0 \u00a0(3)<\/p>\n

We will note Wea\u00a0\u00a0<\/sub>= Je<\/sub>K\u22121<\/sup>Ja<\/sub>T<\/sup>.\u00a0\u03b4e<\/sub>\u00a0represents a vector between actual positions of effector points and positions that effectors have to reach.<\/p>\n

Step III<\/strong>: The final configuration, at the end of\u00a0the time step, is corrected by using the value of the constraint\u00a0response:<\/p>\n

x\u00a0= xfree<\/sup>\u00a0+ K\u22121<\/sup>Ja<\/sub>T<\/sup>\u03bba<\/sub><\/p>\n

Quadratic programming<\/strong><\/h1>\n

In the case of a Quadratic Programming (QP) problem, during the Step II<\/strong>, one sets a QP by minimizing the norm of the vector \u03b4e<\/sub>\u00a0given in (3).<\/p>\n

min( 1\/2 \u03b4T<\/sup>e<\/sub>\u00a0\u03b4e\u00a0<\/sub>) = min( 1\/2 \u0014\u03bba<\/sub>\u0015T<\/sup>\u0014WT<\/sup>ea<\/sub>Wea<\/sub>\u03bba\u00a0<\/sub>+ \u03bba<\/sub>\u0015T<\/sup> \u0014Wea<\/sub>\u03b4e<\/sub>free\u00a0<\/sup>) \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 \u00a0(4)<\/p>\n

subject to min \u2264 \u03bba<\/sub> \u2264 max and \u2212 \u03b5\u00a0\u2264 \u03bbe<\/sub> \u2264 \u03b5<\/p>\n

The size of the QP problem is much smaller than solving the problem in the motion space of equation (2), allowing to solve this problem in real-time. This type of problem can be solved by a variety of methods.\u00a0The matrix of the QP, ie. WT<\/sup>ea<\/sub>Wea,<\/sub> is symmetric. If the number of actuators is equal or less than the size of the effector space, the matrix is also positive-definite. In that case, the solution of the minimization is unique. In the opposite case, i.e when the number of actuators is greater than the degrees of freedom of the effector points, the matrix of the QP is only positive-semidefinite. Consequently, the solution could be non-unique.<\/p>\n

In such case, some QP algorithms are able to find one solution among all possible solutions. A new criterion for the minimization can also be introduced, based on the deformation energy. Indeed, this energy Wdef<\/sub> is linked to the mechanical work of the forces exerted by the actuators. Edef<\/sub> can be computed while evaluating the dot product between \u03bba<\/sub> and the displacements of the actuators \u2206\u03b4a<\/sub> = \u03b4a<\/sub> \u2212 \u03b4a<\/sub>free<\/sup>\u00a0due to the actuator forces Edef<\/sub> = \u03bba<\/sub>T<\/sup> \u2206\u03b4a<\/sub> = \u03bba<\/sub>T<\/sup>Waa<\/sub>\u03bba<\/sub>, with Waa<\/sub>=\u00a0Ja<\/sub>K\u22121<\/sup>Ja<\/sub>T<\/sup>. Yet, matrix Waa<\/sub> is positive-definite if the actuators are placed on different nodes of the FEM or with different directions (i.e\u00a0if there is no linear dependancy between lines of Ja<\/sub>). Thus, we can add this energy in the minimization process by replacing equation (4) by:<\/p>\n

min \u0012(1\/2 \u00a0\u0014\u03bba<\/sub>\u0015T<\/sup>(\u0014<\/sub>WT<\/sup>ea<\/sub>Wea\u00a0<\/sub>+ \u03b5Waa<\/sub>)\u00a0\u03bba\u00a0<\/sub>+ \u03bba<\/sub>\u0015T<\/sup>\u00a0\u0014Wea<\/sub>\u03b4e<\/sub>free<\/sup>) \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 \u00a0(5)<\/p>\n

with \u03b5\u000f chosen sufficiently small so that the deformation energy does not disrupt the quality of the effector positioning. In practice, we use\u00a0\u03b5 \u000f = 1e-3 ||WT<\/sup>ea<\/sub>Wea<\/sub>||\u221e<\/sub> \/ ||Waa<\/sub>||\u221e<\/sub> . Thanks to this modification, the QP matrix becomes positive-definiteand a unique solution of the problem can be found.<\/p>\n

Gauss Seidel<\/strong><\/h1>\n

One can also use a Gauss-Seidel (GS) iterative solver to find a solution among the possible solutions. \u00a0Let\u2019s consider that \u2206\u03b4i<\/sub> provides the shift between the position of the effector and a desired position along direction i<\/em>: \u2206\u03b4i<\/sub> = \u03b4i<\/sub>\u2212\u03b4i<\/sub>desired<\/sup>. We want to find a position on the actuators that deforms the structure so that \u2206\u03b4i<\/sub> = 0.\u00a0During each iteration it<\/em>, the algorithm updates one by one the contribution \u03bbj<\/sub> of each actuator j<\/em>, while freezing the contribution \u03bbk<\/sub> of the other actuators.<\/p>\n

\u03b4j<\/sub> \u2212wjj<\/sub>\u03bbj<\/sub>(it)<\/sup>\u00a0= \u2211<\/span>k=0,k<j-1<\/sub><\/em><\/span> wjk<\/sub>\u03bbk<\/sub>(it)<\/sup>\u00a0+ \u2211<\/span>k=j+1,k<N<\/sub><\/em><\/span>\u00a0wjk<\/sub>\u03bbk<\/sub>(it\u22121)<\/sup>\u00a0+\u03b4j<\/sub>free<\/sup>\u00a0 \u00a0 \u00a0 \u00a0 \u00a0 \u00a0(6)<\/p>\n

As we use a GS algorithm, when doing a local update on actuator j<\/em>, the frozen contribution of the actuators [0 \u2192 j\u22121] comes from the current iteration it<\/em> whereas the contribution\u00a0of actuators [j + 1 \u2192 N] (where N is the total number of actuators) comes from the previous iteration it\u22121<\/em>. The local update of actuator j<\/em> provides a new contribution \u03bbj<\/sub>(it)<\/sup>\u00a0and the effector position \u03b4i<\/sub>, along direction i<\/em>, is easily updated:<\/p>\n

\u03b4i<\/sub> += wij<\/sub>(\u03bbj<\/sub>(it)<\/sup>\u2212\u03bbj<\/sub>(it\u22121)<\/sup>) = wij<\/sub>\u2206\u03bbj<\/sub>\u00a0 \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 \u00a0(7)<\/p>\n

Thus, during the update of j<\/em>, we can measure how the actuator can (or not) reduce the value of \u2206\u03b4i<\/sub>.\u00a0When performing this update in 3 dimensions, direction i<\/em> is alternatively replaced by x, y and z. Imposing a variation of \u03bbj<\/sub> will create a 3D motion of the end effector:<\/p>\n

( \u03b4x \u03b4y \u03b4z )t<\/sup>\u00a0+= ( wxj<\/sub> wyj<\/sub> wzj<\/sub>)t<\/sup>\u00a0* \u2206\u03bbj<\/sub> \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 \u00a0(8)<\/p>\n

We will note\u00a0\u03b4*<\/sub>=\u00a0( \u03b4x \u03b4y \u03b4z )t<\/sup>\u00a0and w*j<\/sub>\u00a0= ( wxj<\/sub>\u00a0wyj<\/sub>\u00a0wzj<\/sub>)t<\/sup>. \u00a0\u00a0The control should try to reduce the shift between actual and desired positions, measured by \u2206\u03b4*<\/sub>. However, the actuator j<\/em> can only move the effector along the direction given by w*j<\/sub>\/||w*j<\/sub>||. Consequently, we search for a value of \u2206\u03bbj<\/sub> so that:<\/p>\n

(w*j<\/sub>\u2206\u03bbj<\/sub> ) \u00b7 ( w*j<\/sub>\/||w*j<\/sub>||\u00a0) = \u2212\u2206\u03b4*<\/sub> \u00b7 ( w*j<\/sub>\/||w*j<\/sub>|| ) \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 (9)<\/p>\n

This value can be obtained by using:<\/p>\n

\u2206\u03bbj<\/sub> = \u2212 (wxj<\/sub>\u2206\u03b4x<\/sub> + wyj<\/sub>\u2206\u03b4y<\/sub> + wzj<\/sub>\u2206\u03b4z<\/sub>) \/||w*j<\/sub>||2<\/sup> \u00a0 \u00a0 \u00a0 \u00a0 \u00a0 (10)<\/p>\n

Given the new value of \u2206\u03bbj<\/sub>, one can update \u03bbj<\/sub>(it)<\/sup>\u00a0and the position \u03b4j<\/sub> of the actuator j<\/em> using the GS equation (6).<\/p>\n

Bibliography<\/strong><\/h1>\n
    \n