Optimizing Parameters for Static Equilibrium of Discrete Elastic Rods with Active-Set Cholesky

To capture bending and twisting of one-dimensional elastic strand structures, Pai [17] introduced the Cosserat theory into graphics, which was later extended to dynamical systems [18]. The Cosserat theory has since been adopted within position-based dynamics [19] and further extended with quaternions [20, 4], volume constraints [21], and projective dynamics [22]. Instead of evolving edge rotations, Shi et al. [23] proposed using a volume-like torsion energy to capture twisting effects. Curvature-based approaches were also presented to capture elastic rod dynamics with a smaller number of degrees of freedom (DOFs) [24, 25]. From the various methods for simulating elastic strands, we selected DER [6, 7] because it can handle general bending and twisting cases with few DOFs.

To achieve sag-free simulations (or equivalently, static equilibrium), one common approach is to enforce zero net forces in the dynamical system. Hadap [26] proposed using inverse dynamics to solve nonlinear force equations for reduced multibody systems, although such methods may not necessarily be applicable to maximal coordinate systems. To solve nonlinear force equations, the Asymptotic Numerical Method (ANM) was employed, achieving faster convergence compared to Newton-type optimizers [27, 28], although it remains unclear how to incorporate box constraints, which help prevent stability problems and violations of physical laws. Curvature-based solvers have also been employed, since they need fewer DOFs to represent elastic rod dynamics. Derouet-Jourdan et al. [29] presented a method for achieving static equilibrium of curvature-based elastic strands [24] and extended their work with frictional contacts [30]. The curvature-based formulations have extensively been used for inverse problems to model, design, and fabricate elastic rods [31, 32, 33, 34].

To achieve static equilibrium in more general contexts, Twigg and Kačić-Alesić proposed minimizing the force norm by optimizing rest shape parameters in mass-spring systems consisting of one-dimensional elements [1]. Takahashi and Batty extended their approach to DER [5], in work closely related to ours. We clarify key differences in Sec II-C. Rest shape optimization has also been applied to two-dimensional sheet-like structures [35] and three-dimensional deformable volumes [36, 37]. To enforce zero-net-force constraints, adjoint methods have been extensively used [2, 38, 12]. While these methods are efficient, their implementation tends to be significantly more complicated because they involve Hessians of objective terms due to differentiation of the nonlinear force constraints. Given that DER Hessian treatments are the subject of ongoing study [12, 23], we formulate our problem without explicitly incorporating the Hessian into the parameter optimization. While finite differences could bypass this complexity [39], they are known to be significantly slower than optimization with analytical gradients.

As an alternative, Hsu et al. [3, 4] proposed a global-local initialization method that first computes global forces to achieve static equilibrium and then adjusts local elements to generate such forces. However, this approach would require excessive stiffening of materials and is applicable only to simulation approaches where local force elements can exclusively be defined. As such, the global-local initialization is not applicable to DER, as discussed by Takahashi and Batty [5].

There are several key differences between our work and the closely related prior work by Takahashi and Batty [5]. Their approach optimizes only the rest shape parameters and may fail to achieve zero net forces. By contrast, our method simultaneously optimizes both the material stiffness and rest shape parameters while enforcing zero net forces as a hard constraint, ensuring static equilibrium at solver convergence. We also reduce the computational cost of parameter optimization in two ways. First, we order the optimization parameters in an interleaved manner to yield a banded system, leveraging the one-dimensional structure of strands, which eliminates the need for matrix reordering in Cholesky solvers. Second, by considering the overlapping force spaces with 4D curvatures [6, 9], we include only two rest curvature variables per vertex, reducing the DOF count in the optimization (see Sec. IV for details). Moreover, instead of relying on penalty-based approaches for box constraints [5], we utilize an active-set-based iterative solver, MPRGP [16], to enforce box constraints precisely. We also propose a new active-set Cholesky preconditioner to accelerate the convergence of MPRGP.

We define the inertia objective and its gradient as follows:

where $\mathbf{q}^{*}=\mathbf{q}^{t}+\Delta t\mathbf{\dot{q}}^{t}+\Delta t^{2}\mathbf{M}^{-1}\mathbf{f}_{\mathrm{ext}}$, $\mathbf{q}^{t}$ and $\mathbf{\dot{q}}^{t}$ denote the generalized positions and velocities at time $t$, respectively, $\Delta t$ time step size, $\mathbf{M}$ a diagonal generalized mass matrix [7, 40], and $\mathbf{f}_{\mathrm{ext}}$ generalized external forces.

We define the stretching energy [7] as $E_{\mathrm{st}}(\mathbf{q})=\sum_{i=1}^{N-2}E_{\mathrm{st},i}(\mathbf{x}_{i},\mathbf{x}_{i+1})$, where $E_{\mathrm{st},i}(\mathbf{x}_{i},\mathbf{x}_{i+1})$ is given by

with $s$ denoting a global constant scalar that scales stiffness parameters, $\bm{\alpha}_{i}$ as the stiffness parameter for stretching on edge $i$, $\mathbf{l}_{i}=\left\lVert\mathbf{x}_{i+1}-\mathbf{x}_{i}\right\rVert_{2}$ representing the length of edge $i$, and $\mathbf{\bar{l}}_{i}$ as its rest length. The gradient can be computed as

where $\mathbf{t}_{i}=\frac{\mathbf{x}_{i+1}-\mathbf{x}_{i}}{\mathbf{l}_{i}}$ denotes the unit tangent vector of edge $i$.

We define the bending objective [6, 9] as $E_{\mathrm{be}}(\mathbf{q})=\sum_{i=1}^{N-2}E_{\mathrm{be},i}(\mathbf{y}_{i})$, where $E_{\mathrm{be},i}(\mathbf{y}_{i})$ is given by

and $\mathbf{y}_{i}=(\mathbf{x}_{i-1}^{T},\bm{\theta}_{i-1},\mathbf{x}_{i}^{T},\bm{\theta}_{i},\mathbf{x}_{i+1}^{T})^{T}\in\mathbb{R}^{11}$, $\bm{\beta}_{i}$ denotes the bending stiffness at vertex $i$, and $\bm{\kappa}_{i}$ and $\bm{\bar{\kappa}}_{i}$ are the four-dimensional curvature and rest curvature, respectively [6, 9]. (While two-dimensional curvatures and rest curvatures can be employed with averaged material frames [7], this method has been shown to inadequately evaluate strand bending due to non-unit averaged material frames [41, 12].) The 11-dimensional gradient of (5) is given by [5]

where $\mathbf{J}_{\mathrm{cu},i}\in\mathbb{R}^{4\times 11}$ denotes the Jacobian of $\bm{\kappa}_{i}$ with respect to $\mathbf{y}_{i}$. Here, $\mathbf{J}_{\mathrm{cu},i}^{T}(\bm{\kappa}_{i}-\bm{\bar{\kappa}}_{i})=\sum_{j=0}^{3}(\bm{\kappa}_{i,j}-\bm{\bar{\kappa}}_{i,j})\nabla\bm{\kappa}_{i,j}$, where $\bm{\kappa}_{i,j}$ denotes the $j$th entry of $\bm{\kappa}_{i}$.

We define the twisting objective [7, 40] as $E_{\mathrm{tw}}(\mathbf{q})=\sum_{i=1}^{N-2}E_{\mathrm{tw},i}(\mathbf{y}_{i})$, where $E_{\mathrm{tw},i}(\mathbf{y}_{i})$ is given by

with $\bm{\gamma}_{i}$ representing the twisting stiffness parameter, and $\mathbf{m}_{i}$ and $\mathbf{\bar{m}}_{i}$ the twist (including reference twist [7, 40]) and rest twist, respectively. The 11-dimensional gradient of (7) is

To handle the nonlinear constraints $\mathbf{c}(\mathbf{p})=0$ in (9), we use ALM [15, 35] and define an augmented Lagrangian objective with a Lagrange multiplier $\bm{\lambda}$ and a penalty parameter $\rho$ as

We then reformulate the parameter optimization (9) as

which we optimize by iteratively solving primal and dual subproblems [15]. The primal subproblem, with fixed $\bm{\lambda}^{k}$ for iteration index $k$ (initialized with $\bm{\lambda}^{0}=0$), is defined as

while the dual subproblem becomes a simple vector update:

Thus, efficiently solving the primal subproblem (12) is crucial for fast parameter optimization.

Given the total force due to the DER objectives, represented as $\mathbf{f}(\mathbf{p})=-\nabla E(\mathbf{q})$, it is natural to define $\mathbf{c}(\mathbf{p})=\mathbf{f}(\mathbf{p})$ to enforce zero net forces. However, in this case, the term $\frac{\rho}{2}\left\lVert\mathbf{c}(\mathbf{p})\right\rVert_{2}^{2}$ in (10) becomes equivalent to a quadratic penalty on the total force, which is numerically ill-conditioned and fails to accurately represent force contributions within the system [5]. Therefore, we define $\mathbf{c}(\mathbf{p})$ as

Consequently, $L(\mathbf{p},\bm{\lambda})$ in (10) can be rewritten as

This objective $\eqref{eq:augmented_lagrangian2}$ is numerically equivalent to the formulation proposed in [5] when considering only the rest shape parameters (excluding stiffness parameters) and setting $\bm{\lambda}=0$. Thus, their approach, which treats $\mathbf{c}(\mathbf{p})$ as a soft constraint, can be regarded as a subset of our method.

We impose box constraints on the optimization parameters to prevent stability problems and to more accurately control the resulting strand dynamics. Specifically, we define $\mathbf{\bar{l}}_{\mathrm{min}}\leq\mathbf{\bar{l}}\leq\mathbf{\bar{l}}_{\mathrm{max}}$ (except for $\mathbf{\bar{l}}_{0}$), $\bm{\bar{\kappa}}_{\mathrm{min}}\leq\bm{\bar{\kappa}}\leq\bm{\bar{\kappa}}_{\mathrm{max}}$, $\mathbf{\bar{m}}_{\mathrm{min}}\leq\mathbf{\bar{m}}\leq\mathbf{\bar{m}}_{\mathrm{max}}$, $\bm{\alpha}_{\mathrm{min}}\leq\bm{\alpha}\leq\bm{\alpha}_{\mathrm{max}}$, $\bm{\beta}_{\mathrm{min}}\leq\bm{\beta}\leq\bm{\beta}_{\mathrm{max}}$, and $\bm{\gamma}_{\mathrm{min}}\leq\bm{\gamma}\leq\bm{\gamma}_{\mathrm{max}}$. Unlike the penalty-based approach [5], which requires safety margins to mitigate the risk of violating physical laws (e.g., negative rest length and material parameters), we handle box constraints precisely using MPRGP [16]. Consequently, we set lower and upper bounds for $\mathbf{\bar{l}}$, $\bm{\alpha}$, $\bm{\beta}$, and $\bm{\gamma}$ as

where $\epsilon$ denotes a small positive value (we set $\epsilon=10^{-10}$). While we set the upper bounds of stiffness parameters to $\infty$ to ensure that static equilibrium is achievable [29, 3], finite values can be used if one prefers to avoid excessively stiff materials that may compromise static equilibrium.

We set the lower/upper bounds for $\bm{\bar{\kappa}}$ and $\mathbf{\bar{m}}$ as

where $\bm{\bar{\kappa}}_{0}$ and $\mathbf{\bar{m}}_{0}$ denote the initial rest curvature and rest twist (before parameter optimization), respectively, $\mu$ and $\eta$ the allowed change for rest curvature and rest twist, respectively, and $\mathbf{e}_{\mathrm{be}}$ and $\mathbf{e}_{\mathrm{tw}}$ vectors of all ones with the same dimensions as $\bm{\bar{\kappa}}$ and $\mathbf{\bar{m}}$, respectively. The stable range of $\bm{\bar{\kappa}}$ and thus $\mu$ can vary depending on simulation settings (e.g., strand geometry and material stiffness) [5]; we use $\mu=1$ by default based on our experiments. Additionally, since changes in rest twist are typically much smaller than those in rest curvature, we set $\eta=\mu/4$ to reduce the number of tunable parameters.

We solve the primal, box-constrained nonlinear minimization subproblem (12) using a Newton-type optimizer. One could approximate the box constraints with a penalty objective to ensure fully unconstrained inner linear systems [15, 5], but penalty methods often require tedious parameter tuning to achieve acceptable results. Even with carefully tuned parameters, these methods may lead to compromises that violate box constraints. Therefore, instead of relying on penalty methods, we incorporate box constraints directly when computing the Newton directions, which is equivalent to solving box-constrained quadratic programs (BCQPs) [42].

The Newton direction $\Delta\mathbf{p}$ at $k+1$ iteration can be computed by solving the following BCQP:

where $\Delta\mathbf{p}_{\mathrm{min}}^{k}$ and $\Delta\mathbf{p}_{\mathrm{max}}^{k}$ denote lower and upper bounds for $\Delta\mathbf{p}$, respectively. We omit $\bm{\lambda}^{k}$ in the argument of $\nabla^{2}L(\mathbf{p})$ as it is independent of $\bm{\lambda}^{k}$ (23). Assuming the ideal step length of one [15], we have $\mathbf{p}_{\mathrm{min}}\leq\mathbf{p}^{k+1}=\mathbf{p}^{k}+\Delta\mathbf{p}^{k+1}\leq\mathbf{p}_{\mathrm{max}}$, giving $\Delta\mathbf{p}_{\mathrm{min}}^{k}=\mathbf{p}_{\mathrm{min}}-\mathbf{p}^{k}$ and $\Delta\mathbf{p}_{\mathrm{max}}^{k}=\mathbf{p}_{\mathrm{max}}-\mathbf{p}^{k}$ [42].

Given the augmented Lagrangian objective (15), we can compute its gradient as

where $\nabla R(\mathbf{p})=\mathbf{W}(\mathbf{p}-\mathbf{p}_{0})$, and we denote the Jacobian of $\mathbf{f}(\mathbf{p})$ with respect to $\mathbf{p}$ by $\mathbf{J}\left(=\frac{\partial\mathbf{f}(\mathbf{p})}{\partial\mathbf{p}}\right)\in\mathbb{R}^{(4N-8)\times(7N-14)}$ (details are provided in Sec. IV-E). Note that we treat $\mathbf{M}$ (which is computed with the initial rest length $\mathbf{\bar{l}}$ [7, 40]) as constant, as $\mathbf{M}$ is introduced to properly scale forces (14).

Because the force formulations are primarily linear with respect to the rest shape and stiffness parameters, the second order derivatives of the forces with respect to $\mathbf{p}$ are mostly zero. An exception is the second order derivatives with respect to $\mathbf{\bar{l}}$; however, these components would be projected to zero to ensure SPD inner systems and valid Newton descent directions [15]. Thus, using the exact Hessian with SPD projections would have only negligible benefits while increasing the implementation complexity and computational cost due to the involvement of third order tensors [5]. Therefore, we ignore the second order derivatives with respect to $\mathbf{\bar{l}}$ and approximate the Hessian in a Gauss-Newton style (excluding $\bm{\lambda}$ from arguments):

As this approximated Hessian is guaranteed to be SPD, we can solve the convex BCQP (20) using MPRGP [16], accelerated with our active-set Cholesky preconditioner (see Sec. V).

Algorithm 1 outlines our parameter optimization method. To accelerate convergence and reduce the risk of converging to undesirable local minima, we use warm starting. We set $\rho=10^{6}$, and use $10^{3}$ and $10^{0}$ for the stiffness and rest shape parts of diagonals in $\mathbf{W}$, respectively, so that we prioritize achieving (in order) zero net forces, small changes in the stiffness parameters, and then rest shape parameters. We perform a single Newton iteration to solve our primal BCQP subproblem (12) as it serves as an inner problem within (11). To guarantee a decrease in the objective, we implement a backtracking line search [15]. The iterations are terminated when $\left\lVert\Delta\mathbf{p}\right\rVert_{2}$ falls below a threshold $\epsilon_{p}\ (=10^{-8})$, iteration count exceeds $k_{\mathrm{max}}\ (=100)$, or backtracking line search fails.

Consider the minimization of a quadratic objective $\frac{1}{2}\mathbf{x}^{T}\mathbf{A}\mathbf{x}-\mathbf{x}^{T}\mathbf{b}$ and corresponding linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$. The associated linear preconditioning system can be expressed as $\mathbf{A}\mathbf{z}=\mathbf{r}$, where $\mathbf{z}$ and $\mathbf{r}$ are vectors corresponding to $\mathbf{x}$ and $\mathbf{b}$, respectively. Adopting an active set $\mathbf{a}$ that indicates whether $\mathbf{x}$ is constrained or not (i.e., with an index $i$, $\mathbf{a}_{i}=1$ or $\mathbf{a}_{i}=-1$ if $\mathbf{x}_{i}$ is constrained by its lower or upper bound, respectively, and $\mathbf{a}_{i}=0$ otherwise) [46], we can define a diagonal selection matrix $\mathbf{S}$ with $\mathbf{S}_{ii}=1-|\mathbf{a}_{i}|$. Given $\mathbf{z}_{u}$ and $\mathbf{z}_{c}$ that denote unconstrained and constrained parts of $\mathbf{z}$, respectively, we need to solve the preconditioning system for $\mathbf{z}_{u}$ while enforcing $\mathbf{z}_{c}=0$ [46]. The preconditioning system can be rewritten (incorporating the necessary constraints [47]) as $(\mathbf{S}\mathbf{A}\mathbf{S}^{T}+\mathbf{I}-\mathbf{S})\mathbf{z}=\mathbf{S}\mathbf{r}$, or given equivalently (by splitting it for $\mathbf{z}_{u}$ and $\mathbf{z}_{c}$) as $\mathbf{S}_{u}\mathbf{A}\mathbf{S}_{u}^{T}\mathbf{z}_{u}=\mathbf{S}_{u}\mathbf{r}$ and $\mathbf{z}_{c}=0$, where $\mathbf{S}_{u}$ denotes a matrix consisting of $\mathbf{S}$’s rows associated with unconstrained variables. One approach to solving these preconditioning systems is to reassemble $(\mathbf{S}\mathbf{A}\mathbf{S}^{T}+\mathbf{I}-\mathbf{S})$ or $\mathbf{S}_{u}\mathbf{A}\mathbf{S}_{u}^{T}$ and factorize it using Cholesky decomposition, but this procedure is costly since active sets can change in each MPRGP iteration [46]. We would rather perform Cholesky decomposition only once and reuse its resulting factor.

Our preconditioning strategy utilizing active sets can also be applied to IC preconditioning. Notably, when the system is banded and fully populated within the band, Cholesky and IC factorization are equivalent, making their preconditioning (even with active sets) identical.

In our parameter optimization, while the system is banded, it is not entirely filled because, e.g., optimization variables for rest curvatures and stretching stiffness are not coupled (i.e., the Hessian (23) lacks off-diagonal elements relating them). Cholesky decomposition is guaranteed to succeed for SPD systems (except for the case where pivot elements become negative due to numerical error [50]), whereas IC factorization may fail unless the system matrix possesses certain properties, such as being an M-matrix [51]. This necessitates reordering the system (i.e., pivoting) potentially breaking the banded structures, adding positive diagonals [51], or reverting to GS [44], ultimately reducing preconditioning effectiveness. Given that approximately 97% of the band is filled, it is acceptable to introduce a small number of fill-ins, considering the IC issues above. We therefore prefer full Cholesky factorization.

We begin by comparing bending formulations using 2D and 4D curvatures to justify our choice of 4D curvatures, given that our parameter optimization employs only two rest curvatures per vertex. The analysis is conducted on a horizontal strand discretized with 30 vertices.

To assess the performance improvement achieved through banded systems and reduced rest curvatures, we experiment with a coil-like strand discretized with 2k vertices, as shown in Figure 4 (we include naive initialization for reference). For a fair comparison with previous work [5], we optimize only the rest shape parameters and exclude box constraints to eliminate the influence of active sets. We compare the following:

1. 1.

   Unbanded system with 4D rest curvatures [5];

2. 2.

   Banded system with 4D rest curvatures;

3. 3.

   Banded system with 2D rest curvatures.

Using the approximate minimum degree (AMD) reordering, the sequential arrangement of the rest shape parameters (rest length, rest curvatures, and rest twist) [5] aligns with the banded system, yielding identical results. The AMD reordering took 2.0 ms while the system solve took 5.0 ms. As the banded system (with 4D rest curvatures) eliminates the need for the reordering, we achieve around $29\%$ performance gain.

Due to the redundant force spaces with rest curvatures, the virtually reduced 2D rest curvatures produce results comparable to those obtained with 4D rest curvatures. The computation times were 2.1ms and 5.0ms with 2D and 4D rest curvatures, respectively. As the DOF count is $4N-8$ and $6N-12$, with non-zero counts per row of $22.0$ and $32.0$ on average, the total number of non-zeros is approximately $88N(=4N\times 22)$ and $192N(=6N\times 32)$ for 2D and 4D rest curvatures, respectively. Thus, the performance ratio $0.42=(2.1/5.0)$ is in good agreement with the ratio of the total number of non-zeros $0.46\approx(88N/192N)$.

To demonstrate the efficacy of our active-set-based BCQP solver, we compare it to a penalty-based one, using a horizontal strand discretized with 30 vertices, as shown in Figure 5 (naive initialization included for reference). In this comparison, we optimize rest shape parameters only and use $\mu=0.4$. For the penalty method, we use a Cholesky solver [5]. As we perform only one Newton iteration to solve each BCQP, if we apply ALM (with Lagrange multipliers initialized to zero) to the BCQP, it becomes equivalent to the penalty method [15].

With the penalty method, it is difficult to accurately control rest curvature changes, allowing rest curvatures to violate the box constraints. Consequently, significant rest shape changes caused the tail end of the strand to flip to the left side after the root was perturbed. In contrast, our active-set method precisely constrains curvature changes retaining the tail end of the strand on the right side. While we may need to compromise static equilibrium, our method better preserves the original shape compared to naive initialization.

The penalty method uses 4 Newton iterations with a total system-solving cost of 0.31 ms (0.078 ms per solve), whereas our method uses 2 Newton iterations with 8.0 MPRGP iterations on average per solve, resulting in a total cost of 0.33 ms (0.17 ms per solve). Although MPRGP needs multiple iterations, our method performs Cholesky factorization (which is costlier than triangular solves) only once per Newton iteration, leading to a moderate cost per MPRGP solve. Furthermore, the strict enforcement of box constraints enhances the convergence of Newton iterations. Consequently, the total cost of our method is comparable to that of the penalty-based approach.

To demonstrate the necessity of optimizing both rest shape and material stiffness parameters and enforcing $\mathbf{c}(\mathbf{p})=0$ as a hard constraint, we compare the following schemes:

1. 1.

   Rest shape only: rest shape only optimization [5];

2. 2.

   Rest shape and material stiffness with penalty: simultaneous optimization of rest shape and material stiffness parameters with enforcement of $\mathbf{c}(\mathbf{p})=0$ via quadratic penalty [15] (which is equivalent to ours with $\bm{\lambda}=0$);

3. 3.

   Rest shape and material stiffness with ALM (ours): rest shape and material stiffness parameter optimization with enforcement of $\mathbf{c}(\mathbf{p})=0$ as a hard constraint via ALM.

We compare our method with various schemes to solve the BCQP subproblem (20) using the scene shown in Figure 5.