Self-Consistent Neural Control: Learning Trajectories for Implicit Equilibrium

Introduction

Physical AI aims to enable robots to perceive, reason, and interact with the physical world. For robots that manipulate compliant objects and structures, many interactions can be modeled through mechanical equilibrium rather than explicit state-transition equations. The system state is therefore obtained by solving an equilibrium problem subject to control-dependent boundary conditions. Efficient learning and control through such implicit models are important for complex robotic interactions and have attracted growing interest [28, 25].

Despite this progress, learning and long-horizon control through implicit equilibrium systems remain challenging, particularly when the equilibrium response varies strongly with the control input. Multiple stable equilibria may co-exist, and the equilibrium branch reached by the system can depend on its deformation history [25]. Deformable-object manipulation provides representative examples. As shown in Figure 1, elastic rods provide a relatively well-behaved, bending-dominated case, although the attained equilibrium branch can still depend on the deformation path. Measuring tapes and thin-walled tubes are more instability-prone: increasing curvature can induce cross-sectional deformation, kinking, and abrupt equilibrium changes through the Brazier effect [21] (Figure 1(B)–(C)). More complex structures, such as Slinky-like systems, exhibit strongly coupled bending, stretching, and large geometric reconfiguration (Figure 1(D)). Although these systems differ mechanically, they illustrate a more general computational challenge: the local sensitivity of an implicit state to control can vary rapidly along the optimization trajectory.

Existing approaches differentiate such systems by unrolling iterative solvers, differentiating explicit learned surrogates, or applying implicit differentiation to converged equilibria or fixed points [24, 2] (Hu et al. 2020). Explicit surrogates provide direct gradients but require representative training data and may accumulate model error over long horizons [24, 17]. Implicit models are commonly evaluated through iterative solvers such as Newton–CG or fixed-point iteration. Differentiating through these iterations requires unrolling and backpropagating through the full solution process [27], resulting in deeply nested computational graphs and substantial computational and memory costs. These limitations motivate methods that differentiate directly through the converged implicit condition without unrolling the internal solver iterations.

Recently, adjoint learning was proposed to retain the implicit forward solve while avoiding costly solver unrolling [25]. Rather than differentiating through the internal iterations of the equilibrium solver, it constructs proxy dynamics whose gradients are determined by the sensitivity of the converged equilibrium state to controllable quantities. To improve computational efficiency over long horizons, adjoint learning further adopts a frozen-tangent approximation, in which higher-order variations of the implicit-state sensitivity along the trajectory are neglected [25].

The accuracy of this approximation, however, relies on the implicit-state sensitivity varying sufficiently slowly along the optimization trajectory. This condition can fail in strongly nonlinear regimes, where small changes in control may produce large variations in the local sensitivity or drive the system toward different solution branches [12] (Mróz and Haftka 1994). The neglected higher-order sensitivity variations can then introduce substantial gradient errors, producing poor search directions and driving optimization toward highly sensitive or undesirable regions of the implicit solution landscape. This exposes a fundamental inconsistency: the gradient approximation assumes locally smooth implicit-state sensitivity, yet the optimization trajectory itself is unconstrained from entering regions where this assumption becomes invalid.

We address this inconsistency with Self-Consistent Neural Control (SCNC), a framework for learning control policies in implicit equilibrium systems with both smooth and rapidly varying equilibrium sensitivities (Figure 2). SCNC is self-consistent in the sense that it regularizes the optimization trajectory to satisfy the local smoothness condition on which its own gradient approximation relies. Rather than assuming that higher-order sensitivity variations are negligible, SCNC explicitly controls their magnitude along the trajectory. We derive a bound relating these variations to the resulting gradient error, making gradient validity a consequence of the optimization process rather than an external assumption. We further develop a practical finite-difference formulation that estimates these sensitivity variations using neighboring converged implicit states to avoid high computational cost.

We validate SCNC across a range of implicit equilibrium systems, including elastic rods, thin-walled tubes, measuring tapes, and Slinky-like structures with neural-network-based energy models. The benchmarks span both smooth bending-dominated responses and strongly instability-prone regimes, with validation in simulation and real robotic manipulation experiments. Compared with adjoint learning and the gradient-free iCEM baseline, SCNC achieves lower task loss and higher optimization success rates while retaining computational efficiency comparable to adjoint learning.

In summary, the main contributions of this work are:

• We characterize theoretically and numerically how rapidly varying implicit-state sensitivities can introduce substantial errors under the frozen-tangent approximation, leading to unreliable gradients along the optimization trajectory.

• We develop Self-Consistent Neural Control (SCNC), derive a bound relating higher-order implicit-state sensitivity variations to gradient error, and introduce a practical finite-difference formulation that estimates and regularizes these variations along the optimization trajectory.

• We validate SCNC across smooth and instability-prone implicit equilibrium systems, including elastic rods, thin-walled tubes, measuring tapes, and Slinky-like structures, demonstrating improved task performance and optimization robustness over adjoint learning and the gradient-free iCEM baseline in both simulation and real-world experiments.

Related Work

Learning and Long-Horizon Control of Implicit Equilibrium Systems

Learning-based control can rely on either explicit dynamics models or implicitly defined system responses. Explicit learned dynamics and differentiable simulators enable gradient-based trajectory optimization and physical rollouts [24] (Hu et al. 2020), but long-horizon optimization can suffer from accumulated model error [17].

For quasi-static systems, control instead induces sequences of nonlinear equilibria. Examples include energy-based shell and rod models [14, 3] and learned fixed-point models such as deep equilibrium networks [2]. With coexisting equilibria, warm-started solves can become path- and branch-dependent [6, 25]. Adjoint learning differentiates such sequences without unrolling the inner solver, but its gradient can become unreliable when local sensitivities vary rapidly [25], which motivates SCNC.

Differentiation and Gradient Reliability Through Implicit Equilibria

Gradients through equilibria can be obtained by unrolling the numerical solver or differentiating the converged optimality condition. Unrolling is general but incurs horizon-dependent memory cost and truncated-solve errors [9, 15]. Implicit differentiation instead solves a linear adjoint system, as used in optimization layers, deep equilibrium networks, and modular frameworks [1, 2, 4], but remains sensitive to Jacobian conditioning and solve accuracy [12]. Error guarantees for approximate implicit or one-step differentiation typically rely on contraction or local-convergence assumptions [13, 5].

In nonlinear mechanics, tangent stiffness can become singular near critical equilibria and continuation may switch branches, requiring specialized sensitivity treatments [20] (Mróz and Haftka 1994). Adjoint learning avoids solver unrolling with a frozen-tangent proxy, but does not control sensitivity variation along instability-prone paths [25]. SCNC instead treats gradient reliability as trajectory-dependent, regularizing the learned trajectory so that sensitivity variation remains controlled and gradient validity becomes an outcome of optimization.

Continuation and Branch Selection near Instabilities

Quasi-static manipulation can be viewed as planning over a manifold of equilibrium configurations [6], while direct parameter continuation fails at simple folds [8]. Arc-length and pseudo-arclength methods track response curves through limit points and snap-through [23, 7, 8]; near bifurcations, critical-point detection and branch-switching methods access secondary branches [11, 26], while deflation reveals additional coexisting solutions [10].

These methods primarily improve equilibrium-branch tracking. SCNC instead optimizes a task-driven control trajectory while regulating the sensitivity variation that governs gradient reliability. Thus, continuation targets forward branch tracking, whereas SCNC targets reliable gradient-based trajectory learning.

Methodology

Adjoint learning method

Consider an implicit equilibrium system whose free equilibrium state is denoted as $x$ and the controllable part is denoted as $z$. $\lambda$ is used to denote the continuation parameter. $x$ and $z$ are functions of $\lambda$. We will use $x(\lambda)$ or $z(\lambda)$ when we try to indicate the function itself, but $x_k$ or $z_k$ for the specific states at $\lambda= k$.

In a normal explicit system, state transitions can be easily expressed in the form of $x_{k+1} = \mathcal{T}(x_k, z_k)$. However, when it comes to implicit systems, we have

Consider the setting in Figure 1(A), $z$ represents the controllable DOF state, which is the state of the rod grabbed by the grippers, $x$ represents the free DOF state, and $\mathcal{G}$ is the implicit equilibrium. The robot should learn how to control $z$ to push $x$ to the target $x^*$ following the rule of $\mathcal{G}$.

In adjoint learning [25], a neural controller parameterized by $\Theta$ takes $\lambda$ as input and generates a sequence of control signals $u_{\Theta}(\lambda)$. The controller drives $z$ through known dynamics $f(z(\lambda), u(\lambda))$, while the corresponding free state $x$ is obtained from the implicit equilibrium condition. Since each equilibrium may require many internal solver iterations, directly backpropagating through all solver steps becomes prohibitively expensive over a long horizon.

Sensitivity is defined as:

where $G_x = \nabla_x\mathcal{G}(x,z)$ and $G_z = \nabla_z\mathcal{G}(x,z)$. In specific dynamics, $\mathcal{G}$ is often defined as $\mathcal{G}(x,z) = \nabla_x E(x,z) - F^{\mathrm{ext}}$ and thus $G_x$ and $G_z$ become tangent stiffness and mixed stiffness respectively. The sensitivity $S(\lambda)$ is never formed explicitly; instead, products $S^{\mathsf{T}}v$ are computed via single adjoint linear solve.

The proxy dynamics described by sensitivity along the homotopy path are:

and the gradient of the loss with respect to the neural network parameters follows from Eq. (3):

where $a(\lambda) = \frac{\partial\mathcal{L}}{\partial x}(\lambda)$ and $g(\lambda) = \frac{\partial\mathcal{L}}{\partial z}(\lambda)$ are proxy system costates. $\mathbf{A}(\lambda) = \frac{\partial f}{\partial z}(z(\lambda), u(\lambda))$ and $\mathbf{B}(\lambda) = \frac{\partial f}{\partial u}(z(\lambda), u(\lambda))$ are two constants for fixed task describing how boundary DOF are updated by control signals. Deduction details of Eq. (3) (4) (5) can be found in technical supplement section 1 or adjoint learning paper [25].

Error from frozen-tangent approximation

Eq. (3) (4) (5) and the transition of costate rely on a frozen-tangent approximation: the equilibrium sensitivity $S(\lambda)$ is treated as locally constant during the backward pass. This amounts to neglecting the variation of $S(\lambda)$ along the homotopy trajectory:

where

and

The frozen-tangent approximation claims that for the system with smoothly varying equilibrium – such as elastic rods with simple quadratic energy – these higher-order variations are often small. Thus these terms can be neglected to avoid calculating third-order derivatives of the energy functional which is really expensive.

However, for instability-prone systems like Slinky-like structures, measuring tapes and thin-walled tubes, both $G_x(\lambda)$ and $G_z(\lambda)$ can vary rapidly along the homotopy path. In these regimes, the neglected higher-order terms in sensitivity can be fatal. Numerical details for detecting $\frac{dS}{d\lambda}$ are included in technical supplement section 7.

Self-consistent learning

Our goal is to constrain the gradient error coming from the frozen-tangent approximation in an acceptable range without any expensive computation.

Note that throughout this paper, $\|\cdot\|$ denotes the Euclidean norm when applied to vectors and the Frobenius norm when applied to matrices or tensors.

Define:

The quantity $\epsilon$ measures the inverse-weighted stiffness variation along the realized homotopy trajectory. We claim that this quantity controls the error introduced by the frozen-tangent approximation up to condition-number constants.

Proof. Define

Based on Eq. (4) (5), the exact and proxy gradients (gradient in adjoint learning method) are

where $i=e$ denotes exact gradient, $i=p$ denotes proxy gradient. $\bar{S}$ denotes the locally frozen sensitivity and $a_e(\lambda)=a_p(\lambda)=a$ (see in technical supplement section 2). Note that $K$ is determined by computation graph of neural network itself and the known dynamics $f$. So $K$ has no difference between exact and proxy gradient. Therefore,

Since uniformly bounded, we can assume $\|I(\lambda)\|\le M_I$, $I\in\{a,K,S\}$. Then

From Eq. (6),

Then

where $\Delta\lambda:=\lambda_{k+1}-\lambda_k$. The transition equation of the costate $\mathbf{g}$ and $\sup_{\lambda\in[0,1]}\|A(\lambda)\|\le M_A$, combined with Gronwall’s inequality give

Details of evaluation on $\mathbf{g}$ are included in technical supplement section 3.

Therefore, Eq. (10) holds with

where $C$ is a constant independent of $\epsilon$.

The premises of theorem 1 are encouraged by the penalty below.

This bound motivates a self-consistent control framework: the controller is optimized not only to minimize the task objective, but also to select trajectories along which the approximation used to compute its gradient remains valid.

Finite-difference surrogate penalty

Theorem 1 suggests constraining the sensitivity variation through

where $\mu_x$ and $\mu_z$ are regularization weights. However, directly evaluating and differentiating Eq. (19) is expensive because of the derivative and inverse matrix based on Eq. (7) (8).

We therefore construct a finite-difference approximation. Discretise the homotopy path into $K$ intervals with $\Delta\lambda=\frac{1}{K}$. The stiffness derivatives can be approximated by

For any matrix $Q$,

where $(\sigma_{\min})$ and $(\sigma_{\max})$ denote the smallest and largest singular values, respectively. Therefore, as long as $G_x$ is not strictly singular, the inverse-weighted quantities in Eq. (22) are controlled. Actually, it is in practical applications and is not only dependent on the configurations but also on Young’s modulus. Across all 24 experimental cases under SCNC, the smallest minimum singular value of $G_x$ observed over all steps is $3.7 \times10^{-4}$, occurring in the thin-walled tube scenario. At this scale, the effect associated with the near-singularity can still be effectively controlled by the $\mu$. We combined the finite-difference method in Eqs. (23)–(24), which motivates the following condition-number-dependent practical surrogate:

The final training objective is

Unlike Eq. (19), Eq. (23) doesn’t contains a matrix inverse. It only penalizes the variation of the tangent and mixed stiffness matrices between consecutive equilibrium states, providing a practical approximation of the original regularization loss.

Pseudocode is included in technical supplement section 8.

Result

We evaluate SCNC using a set of 2D deformable manipulation testbeds spanning smooth and instability-prone implicit equilibrium regimes. The experiments include elastic rods, thin-walled tubes, and measuring tapes in both simulation and the real world, together with simulated Slinky-like structures governed by a learned neural energy model.

We use python jax library to realize our Newton iterative simulator combined with line search when minimizing the energy for implicit equilibrium.

We compare our SCNC with adjoint learning and iCEM. All the methods share the same forward rollout method described above but utilize different ways to calculate gradients.

Experimental setup

We model the elastic rods, thin-walled tubes and measuring tapes as center lines with stretching and bending energy. The center lines are discretized into 101 nodes with each node including 2 DOF: $x$ and $y$. For Slinky-like structures, we follow the Slinky model of (Li et al. 2023) whose energy is governed by a learned neural network. We choose one with 50 coils and each coils has 3 DOF: $x$, $y$ and $\alpha$ which denotes the orientation of the coil. For all of these, the robot arm will grab the both ends. Models’ details are provided in the technical supplement sections 4 and 5.

Tasks and metrics

For all objects, we use the same node-targeting task, in which the robot manipulates both ends of the object to move a designated node to a target position. The control signal $u_{\Theta}(\lambda_i)$

varies from $\mathbb{R}^2$ to $\mathbb{R}^4$, including $x$ axis displacement for two ends and the in-plane rotation for two clamps.

Denote the coordinates of chosen node as $\mathbf{n}$ which is a function of $\lambda$ and target position as $\mathbf{n}^*$. The task loss is defined as:

Besides, we define a optimization is stable if the equilibrium doesn’t switch to the bad branch – kink doesn’t propagate (change the relative position on the object), disappear and no new kink occurs for tapes and tubes; stiffness of Slinky-like structures doesn’t vary too fast (the final distance between left end and right end is less than 1.5 times its original distance).

Each object contains six evaluation cases for every method.

Simulation results

We reported some representative simulation results in Figure 3. The rendered figure shows the stability of our method. As described before, our method is less prone to switching to undesirable equilibrium branches. For case A, adjoint learning gives a bad final equilibrium, where the two initial kinks merge into one and the whole tape twists. While SCNC translates the tape as a whole and preserves the stability of the kink when trying to move the controlled point to the target. For case B, adjoint learning propagates a kink initially near the left end toward the right end of the tape while SCNC keeps the relative position of the kink. Similar things happen on tube cases. For cases C and D, adjoint learning’s trajectories have new kink occurring while SCNC maintains the original configurations and realizes the target at the same time.

When it comes to the Slinky-like structures, SCNC still shows advantages. For case G, adjoint learning primarily moves only the right end of the Slinky, resulting in a substantial increase in its length. This significantly changes the stiffness and sensitivity of the system, making control much more challenging. In contrast, SCNC again learns to translate the entire Slinky toward the target while avoiding excessive stretching. Similar behaviors can be observed in case H.

For elastic rods (cases E and F), the policies learned by adjoint learning and SCNC are nearly indistinguishable, which is consistent with our expectations. Owing to the simple quadratic energy form, elastic rods exhibit smooth and well-behaved sensitivity variations. As a result, the frozen-tangent approximation underlying adjoint learning remains sufficiently accurate, allowing it to achieve performance comparable to that of SCNC.

Baselines comparison

In addition to adjoint learning, we choose iCEM as another baseline because it is widely used for trajectory optimization. [22] We select population size $P = 10$, which already imposes a substantial computational burden. All other parameters needed are included in technical supplement section 6.

Convergence Figure 4(A) shows the score ($\log_{10}$(best loss)) vs. epoch loss history. Each thin curve corresponds to the average loss trajectory of a method for one object, whereas the thick curve represents the overall mean loss trajectory. The results show that both adjoint learning and SCNC substantially outperform iCEM, while SCNC exceeds adjoint learning by 0.601 points in score and achieves an approximately $4\times$ lower average loss at epoch 30.

The largest performance gap between SCNC and adjoint learning is observed for thin-walled tubes, where the average loss differs by more than 1.5 orders of magnitude. Technical Supplement section 7 further shows that this system exhibits the largest sensitivity variations, together with the largest difference in sensitivity variation between the two methods, supporting the effectiveness of explicitly controlling sensitivity variations in instability-prone systems.

In contrast, for elastic rods with relatively small sensitivity variations, adjoint learning already performs well. Consequently, the gap between SCNC and adjoint learning is much smaller, as expected.

Statistical significance analysis is included in technical supplement section 9.

Stability Figure 4(B) compares the stability of each method across different objects according to the stable optimization criterion defined in § 4.2. Among the 24 evaluation cases, SCNC achieves stable optimization in 20 cases (83.3%), substantially outperforming iCEM with 12 cases (50.0%) and adjoint learning with 10 cases (41.7%). In particular, SCNC remains stable on several instability-prone tape and tube tasks where both baselines fail, demonstrating its improved robustness to abrupt equilibrium transitions and rapidly varying sensitivities.

Computational cost As Table 1 reports, SCNC incurs almost no additional overhead compared with adjoint learning, with the average time per epoch increasing by only about 8%. Nevertheless, it remains substantially more efficient than iCEM. This demonstrates the practicality of our finite-difference surrogate, which significantly improves optimization performance while introducing only a marginal computational overhead.

Method	Tape		Tube		Rod		Average
	Time[s]	Mem.[MB]	Time[s]	Mem.[MB]	Time[s]	Mem.[MB]	Time[s]	Mem.[MB]
Adjoint learning	11.56	1834.8	2.85	1865.9	2.00	1680.3	5.47	1793.7
SCNC(ours)	11.77	1849.0	3.26	1873.5	2.76	1762.5	5.93	1828.3
iCEM	113.08	1846.4	21.40	1871.0	17.09	1761.5	50.52	1826.3

Peak memory consumption is comparable across all three methods because the libraries occupy the majority of memory needed. Finally, the computational cost for the Slinky-like structure is reported separately in the Technical Supplement section 6.4 because its energy is represented by a neural network.

Real world experiments

As shown in Figure 5, the real-world experiments follow the simulation setup in Figure 3, except that the Slinky-like structure is excluded and only SCNC results are shown for the elastic rod, where SCNC and adjoint learning perform comparably. Matching subfigure labels denote the same experimental cases.

Consistent with § 4.2.1, the real-world experiments show that SCNC substantially outperforms adjoint learning on thin-walled tubes and measuring tapes, where instability-induced transitions are more pronounced, while achieving comparable performance on elastic rods. This confirms that the advantages of SCNC extend to real-world robotic manipulation.

Additionally, a complete video of each real-world experimental case is provided in the media supplement, using the same label as in Figure 3 and Figure 5.

Discussion

Our results suggest that optimization stability in implicit equilibrium systems is closely tied to the validity of the gradient approximation. When local equilibrium sensitivities vary rapidly, the frozen-tangent approximation can produce inaccurate search directions, leading to unstable optimization and degraded task performance. Rather than recovering the neglected higher-order derivatives, SCNC steers optimization toward trajectories where the approximation remains reliable, resulting in more stable learning and improved performance across instability-prone regimes.

Although demonstrated on energy-based mechanics systems, the proposed framework is not restricted to any particular constitutive model. The derivation only assumes an implicitly defined equilibrium together with differentiable local equilibrium sensitivities, making the idea potentially applicable to a broader class of implicit physical systems. The experiments on neural-network energy models provide an indication of this generality.

The current work nevertheless has several limitations. First, although the proposed theory is generally applicable to implicit equilibrium systems, the experimental evaluation is currently limited to node-targeting manipulation tasks. While these tasks capture the optimization challenges caused by rapidly varying equilibrium sensitivities, they do not fully represent the diversity of practical robotic manipulation problems. Future work will evaluate SCNC on more challenging scenarios, including shape control, 3D manipulation and contact-rich interaction. Second, the regularization weights are currently selected empirically based in the Young’s modulus of the material (See technical supplement section 7). Developing adaptive regularization strategies that automatically adjust the penalty according to the sensitivity variation is another promising direction for future research.

Conclusion

We presented Self-Consistent Neural Control (SCNC), a learning framework for implicit equilibrium systems that improves optimization by explicitly controlling the variation of local equilibrium sensitivity. Unlike adjoint learning, which applies the frozen-tangent approximation with high risk, SCNC regularizes the trajectory to maintain the reliability of the approximation. We further derived a theoretical bound relating sensitivity variation to gradient error and proposed an efficient practical finite-difference surrogate. Experimental results on both physical and learned equilibrium models demonstrate that SCNC achieves more stable optimization and better task performance than existing gradient-based and gradient-free baselines, particularly in instability-prone regimes. We hope this work provides a useful perspective for scalable learning and control in implicit physical systems and inspires future research on robust optimization for physical AI.