Introduction
This section details the formulation that is used for real-time inverse method on FEM. It describes how it can be used to estimate the external loads, pressures, displacements… that lead to a given deformation. The static FEM accounts for non-linear geometrical deformations and integrates over the structure a Hookean constitutive law. During each step i of the simulation, a linearization of the internal forces is computed:
f(xi) ≈ f(xi−1) + K(xi−1)dx (1)
where f provides the volumetric internal stiffness forces at a given position x of 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 − xi−1. 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:
−K(xi−1)dx = p + f(xi−1) + JT λ (2)
where p accounts for the external forces (e.g. gravity) that are known and JT λ gathers the contributions of the Lagrange multipliers. Two types of multipliers are defined:
actuators multipliers λa: we use these constraints to describe the external efforts applied on the boundaries (such as pressure for instance). The location of the boundary conditions are supposed to be known, the directions of the effort JT can be updated at each step, and λa is the unknown intensity of the efforts. We can set (and update at each step i) an interval of potential values as min ≤ λa ≤ max.
effectors multipliers λe: as we don’t want to act on the effectors but only on the actuators, we assign λe = 0. Even if they are null, these multipliers are useful to build the optimization problem.
Indeed, the next step consists in the projection of the FEM model equations into the constraint space: the size of matrix K is often very large so an optimization in the motion space would be very computationally expensive. Instead, using the Schur complement of the constraint problem, we do a projection that drastically reduces the size of the search space. Thus three steps are performed.
Step I: a free configuration xfree of the deformable model is found by solving equation (2) with λ = 0. For the constraints defined on effectors, we compute a violation noted δefree. For instance, in a position effector, δefree provides a vector between the effector position during the free motion and the desired position.
Step II: 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:
δe = JeK−1JaTλa + δefree (3)
We will note Wea = JeK−1JaT. δe represents a vector between actual positions of effector points and positions that effectors have to reach.
Step III: The final configuration, at the end of the time step, is corrected by using the value of the constraint response:
x = xfree + K−1JaTλa
Quadratic programming
In the case of a Quadratic Programming (QP) problem, during the Step II, one sets a QP by minimizing the norm of the vector δe given in (3).
min( 1/2 δTe δe ) = min( 1/2 λaTWTeaWeaλa + λaT Weaδefree ) (4)
subject to min ≤ λa ≤ max and − ε ≤ λe ≤ ε
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. The matrix of the QP, ie. WTeaWea, 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.
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 is linked to the mechanical work of the forces exerted by the actuators. Edef can be computed while evaluating the dot product between λa and the displacements of the actuators ∆δa = δa − δafree due to the actuator forces Edef = λaT ∆δa = λaTWaaλa, with Waa= JaK−1JaT. Yet, matrix Waa is positive-definite if the actuators are placed on different nodes of the FEM or with different directions (i.e if there is no linear dependancy between lines of Ja). Thus, we can add this energy in the minimization process by replacing equation (4) by:
min (1/2 λaT(WTeaWea + εWaa) λa + λaT Weaδefree) (5)
with ε chosen sufficiently small so that the deformation energy does not disrupt the quality of the effector positioning. In practice, we use ε = 1e-3 ||WTeaWea||∞ / ||Waa||∞ . Thanks to this modification, the QP matrix becomes positive-definiteand a unique solution of the problem can be found.
Gauss Seidel
One can also use a Gauss-Seidel (GS) iterative solver to find a solution among the possible solutions. Let’s consider that ∆δi provides the shift between the position of the effector and a desired position along direction i: ∆δi = δi−δidesired. We want to find a position on the actuators that deforms the structure so that ∆δi = 0. During each iteration it, the algorithm updates one by one the contribution λj of each actuator j, while freezing the contribution λk of the other actuators.
δj −wjjλj(it) = ∑k=0,k<j-1 wjkλk(it) + ∑k=j+1,k<N wjkλk(it−1) +δjfree (6)
As we use a GS algorithm, when doing a local update on actuator j, the frozen contribution of the actuators [0 → j−1] comes from the current iteration it whereas the contribution of actuators [j + 1 → N] (where N is the total number of actuators) comes from the previous iteration it−1. The local update of actuator j provides a new contribution λj(it) and the effector position δi, along direction i, is easily updated:
δi += wij(λj(it)−λj(it−1)) = wij∆λj (7)
Thus, during the update of j, we can measure how the actuator can (or not) reduce the value of ∆δi. When performing this update in 3 dimensions, direction i is alternatively replaced by x, y and z. Imposing a variation of λj will create a 3D motion of the end effector:
( δx δy δz )t += ( wxj wyj wzj)t * ∆λj (8)
We will note δ*= ( δx δy δz )t and w*j = ( wxj wyj wzj)t. The control should try to reduce the shift between actual and desired positions, measured by ∆δ*. However, the actuator j can only move the effector along the direction given by w*j/||w*j||. Consequently, we search for a value of ∆λj so that:
(w*j∆λj ) · ( w*j/||w*j|| ) = −∆δ* · ( w*j/||w*j|| ) (9)
This value can be obtained by using:
∆λj = − (wxj∆δx + wyj∆δy + wzj∆δz) /||w*j||2 (10)
Given the new value of ∆λj, one can update λj(it) and the position δj of the actuator j using the GS equation (6).
