**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(x_{i}) ≈ f(x_{i−1}) + K(x_{i−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 = x_{i} − x_{i−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(x_{i−1})dx = p + f(x_{i−1}) + J^{T} λ (2)

where p accounts for the external forces (e.g. gravity) that are known and J^{T} λ 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 J

^{T }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 x^{free} of the deformable model is found by solving equation (2) with λ = 0. For the constraints defined on effectors, we compute a violation noted δ_{e}^{free}. For instance, in a position effector, δ_{e}^{free} 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} = J_{e}K^{−1}J_{a}^{T}λ_{a} + δ_{e}^{free} (3)

We will note W_{ea }= J_{e}K^{−1}J_{a}^{T}. δ_{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 = x^{free} + K^{−1}J_{a}^{T}λ_{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 δ^{T}_{e} δ_{e }) = min( 1/2 λ_{a}^{T}W^{T}_{ea}W_{ea}λ_{a }+ λ_{a}^{T} W_{ea}δ_{e}^{free }) (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. W^{T}_{ea}W_{ea,} 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 W_{def} is linked to the mechanical work of the forces exerted by the actuators. E_{def} can be computed while evaluating the dot product between λ_{a} and the displacements of the actuators ∆δ_{a} = δ_{a} − δ_{a}^{free} due to the actuator forces E_{def} = λ_{a}^{T} ∆δ_{a} = λ_{a}^{T}W_{aa}λ_{a}, with W_{aa}= J_{a}K^{−1}J_{a}^{T}. Yet, matrix W_{aa} 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 J_{a}). Thus, we can add this energy in the minimization process by replacing equation (4) by:

min (1/2 λ_{a}^{T}(_{}W^{T}_{ea}W_{ea }+ εW_{aa}) λ_{a }+ λ_{a}^{T} W_{ea}δ_{e}^{free}) (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 ||W^{T}_{ea}W_{ea}||_{∞} / ||W_{aa}||_{∞} . 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}−δ_{i}^{desired}. 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} −w_{jj}λ_{j}^{(it)} = ∑* _{k=0,k<j-1}* w

_{jk}λ

_{k}

^{(it)}+ ∑

*w*

_{k=j+1,k<N}_{jk}λ

_{k}

^{(it−1)}+δ

_{j}

^{free}(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} += w_{ij}(λ_{j}^{(it)}−λ_{j}^{(it−1)}) = w_{ij}∆λ_{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} += ( w_{xj} w_{yj} w_{zj})^{t} * ∆λ_{j} (8)

We will note δ_{*}= ( δx δy δz )^{t} and w_{*j} = ( w_{xj} w_{yj} w_{zj})^{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} = − (w_{xj}∆δ_{x} + w_{yj}∆δ_{y} + w_{zj}∆δ_{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).