# Closed-forms of planar Kirchhoff elastic rods: application to inverse geometry

Olivier Roussel, Marc Renaud, Michel Taïx

► To cite this version:

Olivier Roussel, Marc Renaud, Michel Taïx. Closed-forms of planar Kirchhoff elastic rods: application to inverse geometry. IMA Conference on Mathematics of Robotics, Sep 2015, Oxford, United Kingdom. 9p. hal-01153456v2

HAL Id: hal-01153456

https://hal.science/hal-01153456v2

Submitted on 20 Feb 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# Closed-forms of planar Kirchhoff elastic rods: application to inverse geometry

By Olivier Roussel, Marc Renaud and Michel Taïx

CNRS, LAAS, Univ. de Toulouse, UPS, 7 avenue du Colonel Roche, F-31400 Toulouse, France

## Abstract

In this paper, we address the problem of inverse geometry for Kirchhoff elastic rods. Based on the explicit formulation of extremal curves in terms of elliptic functions, we derive closed forms for rod shape and sensitivity in the planar case and we show how this can be used efficiently in robotics applications. More specifically, a numerical optimization approach is presented to solve the inverse geometry problem in the general 3-D case and applied to the planar case using closed forms previously introduced.

## 1. Introduction

Consider an elastic rod manipulated by robotic grippers at its both extremities. It is well known that the shape of the rod at equilibrium configurations correspond to minimal energy curves and that there exists a countable number of rod configurations for given grippers placement. Bretl and McCarthy (2014) shown the configuration space for Kirchhoff elastic rods is a 6-dimensional manifold which can be parametrized by one single chart. This result has a lot of potential applications in motion planning such, for example, in assembly and disassembly studies as investigated in Roussel et al. (2015a).

Borum et al. (2014) addressed the problem of state estimation for planar elastic rods, where several markers have to be given to be able to retrieve the corresponding configuration of the rod. This paper can be seen as extension of this work by addressing the problem of finding one of the rod quasi-static configurations that satisfy a given position of only one marker of the rod.

In robotics, the task is usually defined by the end effector of the manipulator. The inverse geometry problem (also named inverse kinematics in the literature) for a kinematic chain such a robotic manipulator arm consists in finding a configuration that satisfy a given end effector position. General solutions for this problem are well known and in some cases closed forms are available (e.g. for 6-R manipulators in Renaud (2006)). For task-based motion planning problems, the inverse geometry must be typically computed more than thousands of times and analytic solutions are critical as they can drastically reduce the computational effort.

When solving the motion planning problem for robotic manipulator arms, deformable cables are usually omitted and the resulting solution path might be inapplicable due to cables constraints. These deformable objects are well studied in the mechanic literature and accurate models are available. However, their integration usually requires finite elements schemes which involve a high computational cost. Assuming quasi-staticity, such cables could be modeled by Kirchhoff elastic rods and the direct and inverse geometry can be computed efficiently, especially in the planar case where, as we will show, closed forms are available to compute rod direct model.

By analogy, this paper focus on the inverse geometry problem for a Kirchhoff elastic rod
by only considering the position of the rod extremity. We investigate a general approach for inverse geometry of 3-dimensional elastic rods based on numerical optimization. To this end, we rely on computing both rod shape and sensitivity that can be done by solving differential systems and implies using costly numerical integration schemes. Alternatively, we give closed forms of rod shape and sensitivity in the planar case and show with experiments this leads to an efficient approach to inverse geometry for planar rods.

Section 2 presents some already known results (see Biggs et al. (2007)) about closed forms of extremal curves for equilibrium configurations, and by treating the planar case as a special case of the 3-dimensional rod, we will derive explicit formulations for rod shape and sensitivity in terms of elliptic functions. In section 3, we investigate approaches based on numerical optimization to solve the inverse geometry problem and apply them with analytic forms in the planar case presented in section 2.

## 2. Explicit formulation of planar elastic rods shape and sensitivity

### 2.1. General case of 3-D rods

Consider an inextensible, non-shearable and unit length linearly elastic rod. The shape of the rod traces a curve that we will describe by the mapping $q: [0, 1] \to SE(3)$. The position along the rod is parametrized by $t \in [0, 1]$ and we will name "base" and "tip" of the rod its extremity at $t=0$ and $t=1$ respectively. Let the mappings $u_1(t), u_2(t), u_3(t)$ such that $u_i: [0, 1] \to \mathbb{R}$ be axial and bending rod strains respectively, and $c_1, c_2, c_3$ be the constants that reflect its elasticity properties. As in Jurdjevic (2005), we say the elastic rod is in static equilibrium in the sense of Kirchhoff if it locally minimizes the elastic energy defined by

$$E_{el} = \frac{1}{2} \int_{0}^{1} \sum_{i=1}^{3} c_i u_i^2 dt.$$

Without loss of generality, we will also assume that the base of the rod is held fixed at the origin, i.e. $q(0) = e$ where $e$ is the identity element of $SE(3)$. Under these assumptions, we will denote by $\mathcal{B}$ the set of positions that the other extremity of the rod $q(1)$ can reach. As shown in Bretl and McCarthy (2014), the problem of static equilibrium of such rods can be formulated as an optimal control problem by

$$\begin{array}{ll}
\underset{q,u}{\text{minimize}} & \frac{1}{2} \int_{0}^{1} \sum_{i=1}^{3} c_i u_i^2 dt \\
\text{subject to} & \dot{q} = q \left( \sum_{i=1}^{3} u_i X_i + X_4 \right) \\
& q(0) = e, \quad q(1) = b
\end{array}
\qquad (2.1)$$

for some $b \in \mathcal{B}$ and where

$$\begin{align*}
X_1 &= \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} & X_2 &= \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} & X_3 &= \begin{bmatrix} 0 & -1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \\
X_4 &= \begin{bmatrix} 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} & X_5 &= \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} & X_6 &= \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \end{bmatrix}
\end{align*}$$

is a basis for $\mathfrak{se}(3)$, the Lie algebra of $SE(3)$. Note that when solving this optimal control problem, the rod tip position $b$ is not an input.

In these conditions, the Maximum Principle states that solutions to this optimal control
EXPLICIT FORMULATION OF PLANAR ELASTIC RODS SHAPE AND SENSITIVITY3

problem are the projections of extremal curves defined on the cotangent bundle $T^*SE(3)$ onto $SE(3)$. Thanks to the Lie Group structure of $SE(3)$, the Hamiltonian can be reduced on the dual of the Lie algebra $sc(3)^*$ and the corresponding (time-varying) Hamiltonian vector fields $\mu: [0, 1] \to sc(3)^*$ can be expressed by

$$ \begin{aligned} \dot{\mu}_1 &= \frac{\mu_3\mu_2}{c_3} - \frac{\mu_2\mu_3}{c_2} \\ \dot{\mu}_2 &= \mu_6 + \frac{\mu_1\mu_3}{c_1} - \frac{\mu_1\mu_3}{c_3} \end{aligned} \qquad \begin{aligned} \dot{\mu}_3 &= -\mu_5 + \frac{\mu_1\mu_2}{c_2} - \frac{\mu_1\mu_2}{c_1} \\ \dot{\mu}_4 &= \frac{\mu_3\mu_5}{c_3} - \frac{\mu_2\mu_6}{c_2} \end{aligned} \qquad \begin{aligned} \dot{\mu}_5 &= \frac{\mu_1\mu_6}{c_1} - \frac{\mu_3\mu_4}{c_3} \\ \dot{\mu}_6 &= \frac{\mu_2\mu_4}{c_2} - \frac{\mu_1\mu_5}{c_1}. \end{aligned} \qquad (2.2) $$

where vector fields $\mu$ are related to controls $u_i$ by $u_i = c_i^{-1}\mu_i$ for $i \in \{1, 2, 3\}$. Let $\mathcal{A}$ be the set homeomorphic to $\mathbb{R}^6$ and $a \in \mathcal{A}$ such that $a_i \triangleq \mu_i(0)$, $i \in \{1, \dots, 6\}$. It has been shown in Bretl and McCarthy (2014) that coordinates in $\mathcal{A}$ offer a global parameterization to the set of static equilibrium configuration for the rod. In other words, $\mathcal{A}$ offers a 6-dimensional configuration space for quasi-static 3-D elastic rods.

Assuming isotropy and normalized elasticity constants such that $c_i = 1$ for $i \in \{1, 2, 3\}$, we have from (2.2) $\dot{\mu}_1 = 0$. Then $\mu_1$ is a constant of motion with $\mu_1 = a_1$ and

$$ \begin{aligned} \dot{\mu}_2 &= \mu_6 \\ \dot{\mu}_3 &= -\mu_5 \end{aligned} \qquad \begin{aligned} \dot{\mu}_4 &= \mu_3\mu_5 - \mu_2\mu_6 \\ \dot{\mu}_5 &= a_1\mu_6 - \mu_3\mu_4 \\ \dot{\mu}_6 &= \mu_2\mu_4 - a_1\mu_5. \end{aligned} \qquad (2.3) $$

The curvature $\kappa$ and the torsion $\tau$ of the curve can be expressed in terms of $\mu$ by

$$ \kappa^2 = \mu_2^2 + \mu_3^2 \qquad \tau = \mu_1 - \frac{\mu_2\mu_5 + \mu_3\mu_6}{\mu_2^2 + \mu_3^2} $$

and, as mentioned in Bretl and McCarthy (2014), (2.3) is equivalent to

$$ 2\ddot{\kappa} + \kappa^3 - 2\kappa(\tau - \lambda_1)^2 = \lambda_2\kappa \qquad \kappa^2(\tau - \lambda_1) = \lambda_3 \qquad (2.4) $$

where the constants of integration are given by

$$ \lambda_1 \triangleq \frac{a_1}{2} \qquad \lambda_2 \triangleq a_2^2 + a_3^2 + 2a_4 - \frac{a_1^2}{2} \qquad \lambda_3 \triangleq \frac{a_1}{2}(a_2^2 + a_3^2) - (a_2a_5 + a_3a_6). $$

Substituting equations in (2.4) and integrating, we obtain

$$ \dot{\kappa}^2 + \frac{1}{4}\kappa^4 + \lambda_3^2\kappa^2 - \frac{\lambda_2}{2}\kappa^2 = \lambda_4 \qquad (2.5) $$

where the constant of integration $\lambda_4$ is given by

$$ \lambda_4 \triangleq a_5^2 + a_6^2 - \frac{1}{4}(a_2^2 + a_3^2)^2 + \frac{1}{2}(a_2^2 + a_3^2)(a_1^2 - 2a_4) - a_1(a_2a_5 + a_3a_6). $$

By making the change of variable $v = \kappa^2$, (2.5) transforms to

$$ \dot{v}^2 + v^3 - 2\lambda_2v^2 - 4\lambda_4v + 4\lambda_3^2 = 0 $$

As stated in Langer and Singer (1984), this equation is in the form $\dot{v}^2 = P(v)$ with $P$ a cubic polynomial. Let $-\alpha_1, \alpha_2, \alpha_3$ be the zeros of the polynomial $P(v)$ such that

$$ -\alpha_1 \le 0 \le \alpha_2 \le \alpha_3. \qquad (2.6) $$

The squared curvature $v$ can be expressed in terms of elliptic functions by

$$ v(t) = \alpha_3 \left( 1 - n \operatorname{sn}^2 \left( rt + \varphi \middle| m \right) \right) \qquad (2.7) $$

where the parameter $m$, the characteristic $n$ and $r$ can be expressed from the polynomial
zeros by

$$m = \frac{\alpha_3 - \alpha_2}{\alpha_3 + \alpha_1} \qquad n = \frac{\alpha_3 - \alpha_2}{\alpha_3} \qquad r = \frac{1}{2} \sqrt{\alpha_3 + \alpha_1}$$

and with the phase $\varphi$ such that $a_2^2 + a_3^2 = \alpha_3(1 - n \operatorname{sn}^2(\varphi|m))$. Note that from (2.6), we have $0 \le m \le n \le 1$.

As outlined in Jurdjevic (2005), it has been shown the Hamiltonian vector fields in (2.3) is integrable and we have proved it can be expressed in the following form

$$\begin{aligned} \mu_2 &= \kappa \sin \psi \\ \mu_3 &= \kappa \cos \psi \end{aligned} \qquad \mu_4 = \frac{1}{2}(\lambda_2 + \frac{a_1^2}{2} - v) \qquad \begin{aligned} \mu_5 &= -\dot{\kappa} \cos \psi + \kappa \dot{\psi} \sin \psi \\ \mu_6 &= \dot{\kappa} \sin \psi + \kappa \dot{\psi} \cos \psi. \end{aligned} \qquad (2.8)$$

where

$$\psi(t) = \lambda_1 t - \frac{\lambda_3}{\alpha_3 r} \left( \Pi(n, \operatorname{am}(rt + \varphi|m)|m) - \Pi(n, \operatorname{am}(\varphi|m)|m) \right) + \psi(0)$$

with $\Pi(n, u|m)$ the elliptic integral of the third kind and $\operatorname{am}(u|m)$ is the Jacobi amplitude.

## 2.2. Planar case

In motion planning, efficient computations of direct and inverse models are critical and closed-forms can remarkably improve the computational effort. Although neither the curve $q(t)$ nor the rod sensitivity $\frac{\partial q(t)}{\partial a}$ can be explicitly expressed in the general 3-D case, we will show in this section that closed forms can be obtained in the planar case.

### 2.2.1. Curvature and internal wrenches

Considering only planar curves $q(t)$ in the $xy$-plane with $q = (0, 0, \theta, x, y, 0)^T$, Hamiltonian vector fields defined in (2.2) simplify to

$$\begin{aligned} \dot{\mu}_1 &= 0 \\ \dot{\mu}_2 &= 0 \end{aligned} \qquad \begin{aligned} \mu_3 &= -\mu_5 \\ \mu_4 &= \mu_3\mu_5 \end{aligned} \qquad \begin{aligned} \dot{\mu}_5 &= -\mu_3\mu_4 \\ \dot{\mu}_6 &= 0. \end{aligned} \qquad (2.9)$$

Closed-forms of rod internal wrenches $\mu(t)$ defined in (2.8) reduce to

$$\begin{aligned} \mu_1 &= 0 \\ \mu_2 &= 0 \end{aligned} \qquad \begin{aligned} \mu_3 &= \kappa \\ \mu_4 &= -\frac{1}{2}(\kappa^2 + \lambda_2) \end{aligned} \qquad \begin{aligned} \mu_5 &= -\dot{\kappa} \\ \mu_6 &= 0. \end{aligned} \qquad (2.10)$$

And constants of integration defined in (2.1) simplify to

$$\lambda_1 = 0 \qquad \lambda_2 = a_3^2 + 2a_4 \qquad \lambda_3 = 0 \qquad \lambda_4 = a_5^2 - a_3^2(\frac{1}{4}a_3^2 + a_4).$$

We retrieve the same results as we would have obtained by applying the same problem formulation on the Lie Group $SE(2)$ rather than $SE(3)$. Therefore, in the rest of this section we will restrict to solutions of (2.1) that are similar to trajectories on $SE(2)$, which are generated by the subset of initial conditions $\{a \in \mathcal{A} : (a_1, a_2, a_6) = (0, 0, 0)\}$.

In the following equations, when referring to an elliptic function $pq$, we will simplify the notation $pq(u|m)$ to $pq u$. Also, let $\Gamma(t) \triangleq rt + \varphi$, $\Gamma_0 \triangleq \Gamma(0)$ and $\varepsilon \triangleq \operatorname{sgn}(a_3)$.

We can distinguish three cases as outlined in Langer and Singer (1984) and Singer (2007):

* Case I: $\lambda_4 > 0$

$$\kappa(t) = \varepsilon \sqrt{\alpha_3} \operatorname{cn} \Gamma(t)$$

The curvature $\kappa(t)$ oscillates between $\sqrt{\alpha_3}$ and $-\sqrt{\alpha_3}$ and the resulting curve $q(t)$ is called a "wavelike" elastica.
## EXPLICIT FORMULATION OF PLANAR ELASTIC RODS SHAPE AND SENSITIVITY

* Case II: $\lambda_4 < 0$

$$ \kappa(t) = \varepsilon\sqrt{\alpha_3} \operatorname{dn} \Gamma(t) $$

The curvature $\kappa(t)$ is non-vanishing and the curve $q(t)$ is called a "orbit-like" elastica.

* Case III: $\lambda_4 = 0$

$$ \kappa(t) = \varepsilon\sqrt{\alpha_3} \operatorname{sech} \Gamma(t) $$

This corresponds to the borderline case where the curvature is non-periodic.

Note that curvature can be reduced to a unique formulation by assuming $m \in [0, \infty)$ and applying the Jacobi's real transformation (see Abramowitz and Stegun (1964) §16.11).

### 2.2.2. Rod shape $q(t)$

From the differential system defined in (2.1), it follows that

$$ \dot{\theta} = u_3 \qquad \dot{x} = \cos \theta \qquad \dot{y} = \cos \theta. $$

Using (2.10), the integration of the curvature is given by

$$ \begin{gathered} \cos \theta(t) = \beta_1(0)\beta_1(t) + 4\beta_2(0)\beta_2(t) \qquad \sin \theta(t) = 2\varepsilon(\beta_1(0)\beta_2(t) - \beta_2(0)\beta_1(t)) \\ x(t) = \beta_1(0) \int \beta_1(t) + 4\beta_2(0) \int \beta_2(t) \qquad (2.11) \\ y(t) = 2\varepsilon \left( \beta_1(0) \int \beta_2(t) - \beta_2(0) \int \beta_1(t) \right). \end{gathered} $$

The functions $\beta_1(t)$ and $\beta_2(t)$ can be explicitly given using Jacobi elliptic functions and the elliptic integral of second kind $E(u|m|$ in the three cases as follows

* Case I: $\lambda_4 > 0$

$$ \beta_1(t) \triangleq 2 \operatorname{dn}^2 \Gamma(t) - 1 \qquad \beta_2(t) \triangleq \sqrt{m} \operatorname{sn} \Gamma(t) \operatorname{dn} \Gamma(t) $$

$$ \int \beta_1(t) = \frac{2}{r} (\mathrm{E}(\mathrm{am} \Gamma(t)) - \mathrm{E}(\mathrm{am} \Gamma_0)) - t \qquad \int \beta_2(t) = -\frac{1}{r} (\mathrm{cn} \Gamma(t) - \mathrm{cn} \Gamma_0) $$

* Case II: $\lambda_4 < 0$

$$ \beta_1(t) \triangleq 1 - 2 \operatorname{sn}^2 \Gamma(t) \qquad \beta_2(t) \triangleq \operatorname{sn} \Gamma(t) \operatorname{cn} \Gamma(t) $$

$$ \int \beta_1(t) = \frac{1}{m} \left( t(m-2) + \frac{2}{r} \left( \mathrm{E}(\mathrm{am} \Gamma(t)) - \mathrm{E}(\mathrm{am} \Gamma_0) \right) \right) \qquad \int \beta_2(t) = -\frac{1}{r} \left( \mathrm{cn} \Gamma(t) - \mathrm{cn} \Gamma_0 \right) $$

* Case III: $\lambda_4 = 0$

$$ \beta_1(t) \triangleq 2 \operatorname{sech}^2 \Gamma(t) - 1 \qquad \beta_2(t) \triangleq \operatorname{sech} \Gamma(t) \operatorname{tanh} \Gamma(t) $$

$$ \int \beta_1(t) = \frac{2}{r} \left( \tanh \Gamma(t) - \tanh \Gamma_0 \right) - t \qquad \int \beta_2(t) = -\frac{1}{r} \left( \operatorname{sech} \Gamma(t) - \operatorname{sech} \Gamma_0 \right). $$

### 2.2.3. Rod sensitivity

The elastic rod sensitivity is given by the 6-dimensional Jacobian matrix $J(t, a)$ with $[J]_{ij} = \frac{\partial q_i}{\partial a_j}$ in the 3-D case. As we consider only the planar case in this section, we will only focus on the 3-dimensional block of $J(t, a)$ for $i, j \in \{3, 4, 5\}$, i.e. $\frac{\partial \theta}{\partial a_j}$, $\frac{\partial x}{\partial a_j}$ and $\frac{\partial y}{\partial a_j}$ for $j \in \{3, 4, 5\}$. These can be obtained by differentiation of the explicit formulation of the curve $q(t)$ given in (2.11) and using that

$$ \frac{\partial \theta}{\partial a} = \frac{\partial \sin \theta}{\partial a} \cos \theta - \frac{\partial \cos \theta}{\partial a} \sin \theta. $$
The complete differentiation of the curve in terms of elliptic functions can be found in Roussel et al. (2015b).

### 3. Inverse geometry

Although closed forms of rod shape can be obtained in the planar case only, there is no explicit formulation of the inverse geometry problem considered in this paper. Alternatively, we will present numerical approaches to tackle this problem.

#### 3.1. A general numerical approach to inverse geometry

Sufficient conditions for optimality state that the curve $q(t)$ obtained from initial conditions $a$ correspond to a stable equilibrium of the elastic rod if and only if the Jacobian matrix $J(t, a)$ is non-singular for any $t \in (0, 1]$. We will denote the set of initial conditions that satisfy this property by $A_{stable}$. Let $\Upsilon: A \to B$ be the map that transforms initial conditions to the rod tip position $q(1)$ and which Jacobian matrix is $J(1, a)$. If we denote $\Upsilon|_{stable}$ the restriction of $\Upsilon$ by

$$ \Upsilon|_{stable}: A_{stable} \to B_{stable} $$

then by definition $J(1, a)$ is non-singular and $\Upsilon|_{stable}$ is a local diffeomorphism.

We can define now the general problem of inverse geometry for an elastic rod. Given a position of the rod tip $b^{des} \in B_{stable}$, the problem of inverse geometry for an elastic rod consists in finding one of the parameterization $a^* \in A_{stable}$ that satisfies $\Upsilon|_{stable}(a^*) = b^{des}$. Let now $f$ by $f(a) \triangleq \Upsilon|_{stable}(a) - b^{des}$. The presented inverse geometry problem is then equivalent to finding one of the zeros of the function $f$. This can be solved using numerical optimization approaches such as Newton's method as described in Nocedal and Wright (2006). The principle consists in iteratively approximate a model $M_k(p_k)$ of $f(a_k + p_k)$ where $p_k \in A$ is the descent step at the iteration $k$ and to deduce the step $p_k$ for which $M_k(p_k) = 0$. As $f$ is also a local diffeomorphism, we can approximate it linearly using Taylor's theorem with

$$ f(a + p) = f(a) + J(1, a)p + O(\|p\|^2) \quad (3.1) $$

Using this model to compute the Newton's step leads to $p_k = -\alpha J(1, a_k)^{-1}f(a_k)$ for any non-singular Jocabian $J(1, a_k)$, i.e. for $a_k \in A_{stable}$. The step length parameter $\alpha$ is typically chosen in $(0, 1]$. Then, the solution $a^*$ can be iteratively approximated by $a_{k+1} = a_k + p_k$ for a given initial guess $a_0 \in A_{stable}$.

#### 3.2. Experiments in the planar case using closed-forms

The numerical optimization scheme requires several evaluations of the function $\Upsilon(a)$ and the Jacobian $J(1, a)$. In Bretl and McCarthy (2014), solutions are expressed as non-linear differential systems which can be solved through costly numerical integration. On the other hand, in the case of planar rods, we derived in section 2.2 closed forms of $\Upsilon(a)$ (equations (2.11) for $t = 1$) and $J(1, a)$. We implemented both numerical and analytic C++ approaches to compute these maps. The numerical integration was done with a 4-th order Runge-Kutta scheme and requires a sufficiently small time step to minimize the integration error (time step of $10^{-4}$ for a relative error of $10^{-3}$). Consequently, in the planar case $\Upsilon(a)$ and $J(1, a)$ can be computed $10^2$ to $10^3$ faster with closed-forms. Steps of the solving process based for the same rod tip placement $b^{des}$ with three different initial guesses $a_0$ are illustrated in figure 1.

Experimental results for planar rods inverse geometry using analytical forms with
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Planar rod inverse geometry using Newton's method. The same rod tip position $b^{des}$ is used with three different initial guesses $a_0$, leading to different solutions $a^*$ corresponding to "wavelike" configuration(top) or "orbit-like" configurations (middle and bottom).</figcaption></figure>

<table><caption>TABLE 1. Results for inverse geometry in the planar case using analytic forms and Newton's method for various values of the length step α</caption><thead><tr><th>α</th><th>Success (%)</th><th>Unstability failure (%)</th><th>Resolution time (μs)</th><th>Number of iterations</th></tr></thead><tbody><tr><td>1.0</td><td>81.5 ± 15.9</td><td>18.4 ± 15.9</td><td>71.2 ± 7.4</td><td>8.1 ± 0.8</td></tr><tr><td>0.9</td><td>86.9 ± 15.1</td><td>13 ± 15.1</td><td>92.2 ± 7</td><td>10.6 ± 0.8</td></tr><tr><td>0.7</td><td>92.8 ± 13.4</td><td>7.2 ± 13.3</td><td>130.7 ± 8.6</td><td>14.9 ± 0.8</td></tr><tr><td>0.5</td><td>95.7 ± 11.6</td><td>4.2 ± 11.5</td><td>197.9 ± 9.7</td><td>22.8 ± 1.1</td></tr><tr><td>0.3</td><td>97.3 ± 9.98</td><td>2.5 ± 9.9</td><td>361.7 ± 14.3</td><td>41.6 ± 1.6</td></tr><tr><td>0.2</td><td>97.4 ± 9.7</td><td>2.1 ± 9.3</td><td>564.8 ± 15.8</td><td>64.9 ± 1.9</td></tr></tbody></table>

Newton's method are given in table 1. The solver was run on 1000 different values of $b^{bes}$, with 100 different values of the initial guess $a_0$ and with various step length $\alpha$. The solver could fail either by reaching the maximum number of iterations (here fixed at 100) or by falling on an instability where the Jacobian matrix is singular. The stopping criterion was set as an angular and position error of 1.e-6 on the desired tip position $b^{des}$.

These results reflect interesting aspects of the problem. First, a relatively low number of iterations are required for high values of $\alpha$ with a good success ratio. This suggests a low error in $O(\|p\|^2)$ in (3.1) so the choice of a linear approximation for the model is sufficient. Consequently, using a line search to determine the appropriate value of $\alpha$ at
Closed-forms of planar Kirchhoff elastic rods: application to inverse geometry 8

each iteration could slightly reduce the total number of iterations but the computational
cost induced makes this approach inefficient.

As α decreases, the success ratio tends to 100%. This is due to the fact that smaller
steps performs better with non-linearities making the chance to fall on instability de-
creases, but at the cost of increasing the number of iterations. Reciprocally, we could
interpret that for high values of the step length α, the probability to reach the instability
region is getting higher as the problem tends to be highly non-linear when getting closer
to instability. This phenomena can be observed in figure 3.2 (bottom case) at the first
iterations and obviously depends on the stability of the initial guess a₀. It has been shown
in Borum and Bretl (2015) that the space of stable initial conditions Astable is simply
connected and quasi star-shaped under a scaling transformation. This indicates that it
is always possible to find any solution a* for a given position bdes and from any initial
guess a₀ avoiding degenerate cases by staying in Astable. In this direction, the descent
step pₖ should be adjusted to maximize the descent and getting as far as possible from
the instability region. From a general perspective, despite its simplicity the Newton's
method offers sufficient results to be exploited in motion planning applications.

**4. Conclusion and perspectives**

† In this paper, we have presented a general approach to address the problem of inverse geometry for Kirchhoff elastic rods. By taking advantage of analytic forms of rod shape and sensitivity in the planar case, we were able to solve efficiently this problem using numerical optimization techniques. Although closed forms cannot be obtained for rod shape and sensitivity in the 3-D case, we plan to use similar optimization schemes by using both closed forms of extremal curves and numerically approximated shape and sensitivity to obtain reasonable performances. Furthermore, even if preliminary experimental results in the planar case are encouraging, a deeper investigation could improve the general efficiency by considering the distance to instability. Also, we would like to extend the problem addressed by promoting solutions that minimizes the total elastic energy.

ACKNOWLEDGMENT

This work was partially support by the French National Research Agency under the
project Flecto (ANR- Digital Models) and by ERC Grant 340050 Actanthrope.

REFERENCES

M. Abramowitz and I. Stegun. *Handbook of Mathematical Functions*. Dover, 1964.

J. Biggs, W. Holderbaum, and V. Jurdjevic. Singularities of optimal control problems on some 6-d lie groups. *Automatic Control, IEEE Transactions on*, 52(6):1027–1038, June 2007.

A. Borum and T. Bretl. The free configuration space of a kirchhoff elastic rod is path-connected. In *2015 IEEE Int. Conf. on Robotics and Automation (to appear in)*, ICRA 2015, Seattle, USA, May 26 - 30, 2015, 2015.

A. Borum, D. Matthews, and T. Bretl. State estimation and tracking of deforming planar elastic rods. In *2014 IEEE Int. Conf. on Robotics and Automation, ICRA 2014*, Hong Kong, China, May 31 - June 7, 2014, pages 4127–4132, 2014.

T. Bretl and Z. McCarthy. Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations. *I. J. Robotic Res.*, 33(1):48–68, 2014.

V. Jurdjevic. *Integrable Hamiltonian systems on complex Lie groups*. Memoirs of the American Mathematical Society, 2005.
J. Langer and D. A. Singer. The total squared curvature of closed curves. *J. Differential Geom.*, 20(1):1-22, 1984.

J. Nocedal and S. Wright. *Numerical optimization*. Springer Series in Operations Research and Financial Engineering. 2006.

M. Renaud. Calcul des modles gomtriques inverses des robots manipulateurs 6r. Technical Report 06332, LAAS-CNRS, 2006.

O. Roussel, A. Borum, M. Taix, and T. Bretl. Manipulation planning with contacts for an extensible elastic rod by sampling on the submanifold of static equilibrium configurations. In *2015 IEEE Int. Conf. on Robotics and Automation (to appear in), ICRA 2015, Seattle, USA, May 26 - 30, 2015*, 2015a.

O. Roussel, M. Renaud, and M. Taix. Closed-forms of kirchhoff elastic rods shape and sensitivity in the planar case. Technical Report 15083, LAAS-CNRS, 2015b.

D. A. Singer. Lectures on elastic curves and rods, 2007.