Dynamic trajectory tracking for robots involves (i) defining the practical workspace of the robot, (ii) designing a target trajectory in this workspace, and (iii) finding optimal control inputs that command the robot to move along the prescribed trajectory. Tasks (i) and (ii) are specific to the robots and the operating domain of interest, whereas task (iii) involves a generally applicable method. To this end, model predictive control (MPC) is applied over a user-defined planning horizon, taking (i) and (ii) as inputs to the optimization.

We fix the workspace as \( \mathbb{R}^{w}\) with some positive integer \( w\) . We then define \( \mathbf{z} \in \mathbb{R}^{w}\) as the coordinates of the workspace and \( \boldsymbol{\Gamma}(t) \subset \mathbb{R}^{w}\) as the target trajectory defined on a finite time interval \( \mathcal{T}=[t_i,t_f]\) . This interval is divided into smaller, equally spaced planning horizons of length \( \Delta t\) . We therefore have \( \mathcal{T} = \bigcup_{j=0}^{N-1}(t_j,t_{j+1}]\) where \( t_0 = t_i\) , \( t_{N} =t_f\) , \( t_{j+1}-t_{j}=\Delta t\) , and \( N = \frac{t_f -t_i}{\Delta t}\) . In each planning interval \( (t_j,t_{j+1}]\) , we have to solve the following optimization problem :

where \( \mathbf{x} \in \mathbb{R}^n\) is the phase space of the positions and velocities of the robot, \( \mathbf{y} \in \mathbb{R}^p\) is the observable space of the robot and \( \mathbf{u} \in \mathbb{R}^{n_u}\) are the control inputs of the robot. \( \mathcal{Z}\) is a set of constraints in the workspace, and \( \mathcal{U}\) is a compact set of admissible control inputs. \( \mathbf{Q_z}\in \mathbb{R}^{w \times w}\) is a positive semi-definite matrix that represents the workspace variables’ cost and \( \mathbf{R_u}\in \mathbb{R}^{n_u \times n_u}\) is a positive definite matrix that penalizes the control input. The function \( \mathbf{h}:\mathbb{R}^{n} \to \mathbb{R}^{p}\) maps the phase space to the observable space. The function \( \mathbf{c}:\mathbb{R}^{p} \to \mathbb{R}^n\) maps the observable space to the phase space and is only guaranteed to exist and be invertible when \( p=n\) . Lastly, the operator \( \mathbf{C}: \mathbb{R}^{p} \to \mathbb{R}^{w}\) maps the observable space to the workspace.

The initial condition \( \mathbf{y}(t_j)\) is the observed configuration of the robot, inferred from a camera at a time instant \( t_j\) , after applying the optimal control inputs that solve the optimal control problem in the previous planning horizon \( (t_{j-1},t_j]\) . For the first planning horizon \( (t_0,t_1]\) , we define the starting configuration of the robot as \( \mathbf{y}(t_0) = \mathbf{y}_0\) . The quantity \( \mathbf{y}(t_j)\) acts effectively as a feedback to the optimal control problem defined on the target’s horizon \( \mathcal{T}\) . Hence, an optimal control input \( \mathbf{u}(t)\) for the horizon \( \mathcal{T}\) is a solution to a closed-loop control problem. If one decides not to use the observed configurations and rely only on the trajectory of the model \( \mathbf{x}(t)\) , then the optimization problem becomes an open-loop problem.

In practical applications, the robot’s observability is limited, and its actual dynamics are unknown. Models for those dynamics are at best complicated partial differential equations (PDEs). One numerically approximates such a PDE with a very large finite-dimensional set of ordinary differential equations (ODE). Optimizing such a high-dimensional ODE over a given planning horizon is computationally expensive, necessitating a reduced-order model. In the following sections, we will rigorously construct such reduced-order models for soft robots in practical operating domains of interest. We list key symbols and some notation that appear in this paper in Table 1.

General robot motion is governed by an ODE

where \( \mathbb{R}^n\) is the phase space of the robot, \( \mathbf{u}(t) \in \mathbb{R}^{n_u}\) is the control input, and \( \mathbf{f}\) and \( \mathbf{g}\) are nonlinear functions. Setting \( \mathbf{u}(t)=0\) in eq.(2), results in an autonomous dynamical system describing the motion of an uncontrolled robot:

By design, most physical robots without control have at least one equilibrium \( \mathbf{x}^*\) defined by \( \mathbf{f}(\mathbf{x}^*)=0\) . We can always perform a coordinate transformation to shift this equilibrium to the origin, and hence we will assume \( \mathbf{x}^{*}=0\) . We also assume this equilibrium to be asymptotically stable, which translates to \( \text{Re}[\text{Spec}(\mathbf{D}_{\mathbf{x}} \mathbf{f}(0))]<0\) for the spectrum of \( \mathbf{D}_{\mathbf{x}} \mathbf{f}(0)\) . These assumptions hold in our setting, since soft robots, by construction, have an asymptotically stable rest position.

Under appropriate nonresonance conditions among the eigenvalues of \( \mathbf{D}_{\mathbf{x}} \mathbf{f}(0)\) , a hierarchy of spectral submanifolds (SSMs) anchored to \( \mathbf{x}^* =0\) exists (see [28]). In our setting here, (primary) SSMs are the unique smoothest invariant manifolds of system (3) that are tangent to spectral subspaces of \( \mathbf{D}_{\mathbf{x}} \mathbf{f}(0)\) at the origin and have the same dimensions as those subspaces. Secondary SSMs with the same invariance and tangency properties also exist, but they are strictly less smooth than the primary SSMs. In this paper, we will only rely on the primary SSMs and refer to them simply as SSMs for brevity. Specifically, slow SSMs are tangent to spectral subspaces spanned by a set of slowest decaying eigenmodes. A unique slow SSM captures most of the dominant dynamics of the system and hence is an ideal candidate for model reduction in system (3). A key feature that makes these structures robust for control or forced response prediction is that they are normally attracting, which means that any normal rates of attraction to the SSM overpower all the tangential rates of attraction within the SSM.

When the control input \( \mathbf{u}(t)\) is non-zero and uniformly bounded, generalized time-varying SSMs anchored to generalized stationary states of the same stability type as the uncontrolled origin are proven to exist (see [30]). Those generalized SSMs can be rigorously constructed for the system (2) for moderate control inputs. They may not survive for large control inputs, but one can still devise systematic control strategies based on them.

To address this limitation on the control input magnitude, [30] also discuss an alternative case wherein time-dependent SSMs called adiabatic SSMs (aSSMs), exist under slowly varying inputs \( \mathbf{u}(t) = \mathbf{u}^s(\epsilon t)\) for \( \epsilon \ll 1\) in system (2). This adiabatic setting allows for larger control inputs, but assumes that control inputs vary in time at a rate that is much slower than the internal decay rates of the uncontrolled robot to its fixed point. Using the formulation of [30], eq.(2) subjected to those slow control inputs can be re-written as the slow-fast system

Setting \( \epsilon =0\) in these, we obtain

This is an autonomous system of ODEs similar to (3) but describes the dynamics of the robot subjected to constant control inputs \( \mathbf{u}(t) \equiv \mathbf{\bar{u}}\) . We will assume that system (6) admits a critical manifold \( \mathcal{L}_0\) (i.e; a manifold of fixed points) of dimension \( n_u\) in the form

We further assume that for any fixed \( \mathbf{u}(t) \equiv \mathbf{\bar{u}}\) , the fixed point \( \mathbf{\bar{x}} = \mathbf{S}(\mathbf{\bar{u}})\) of the first equation in (4) is asymptotically stable. Under this assumption, we let \( \mathbf{A}(\mathbf{u}) = \mathbf{D}_{\mathbf{x}} \left[\mathbf{f}(\mathbf{x}) + \mathbf{g}(\mathbf{x},\mathbf{u})\right]|_{\mathbf{x}=\mathbf{S}(\mathbf{u})}\) and \( \text{Spec}(\mathbf{A}(\mathbf{u})) = {\lambda_j(\mathbf{u})}_{j=1}^{n}\) . We further assume that for some non-zero negative real number \( K\) , strict spectral splitting exists:

The first \( d\) eigenvalues correspond to slow eigenvectors that span a slow spectral subspace denoted \( E(\mathbf{u})\) . From [28], a unique smoothest \( d\) -dimensional SSM, \( \mathcal{W}(E(\mathbf{u}))\) , exists tangent to \( E(\mathbf{u})\) and anchored to \( \mathbf{x} = \mathbf{S}(\mathbf{u})\) , provided that appropriate nonresonance conditions hold among the eigenvalues \( \lambda_j(\mathbf{u})\) for all \( \mathbf{u} \in \mathcal{U}\) . Note that the collection of these SSMs, \( \mathcal{A}_0 = \cup_{\mathbf{u} \in \mathcal{U}} \mathcal{W}(E(\mathbf{u}))\) , is a \( \rho\) -normally attracting compact invariant manifold, with

as defined by [39].

Pan/Zoom

Pan/Zoom

Pan/Zoom

For \( \epsilon >0\) , by Theorem 6 of [30], the manifold \( \mathcal{A}_0\) survives as a nearly diffeomorphic invariant manifold \( \mathcal{A}_{\epsilon}\) , which we call an adiabatic SSM (aSSM). The critical manifold \( \mathcal{L}_0\) will also survive as a slow trajectory \( \mathcal{L}_{\epsilon}\) that acts as an anchor to the aSSM. We can seek \( \mathcal{L}_{\epsilon}\) and \( \mathcal{A}_{\epsilon}\) as Taylor expansions in \( \epsilon\) using the governing equations. Let \( \mathbf{V}(\mathbf{u}) \in \mathbb{R}^n \times \mathbb{R}^{d}\) be a matrix whose columns form the basis of the spectral subspace \( E(\mathbf{u})\) , and let

denote a local coordinate along \( E(\mathbf{u})\) . Note that \( \mathcal{A}_0\) can be parametrized as \( \mathcal{A}_0 = {\mathbf{x}\in \mathbb{R}^n: \mathbf{x} = \mathbf{W}(\mathbf{r},\mathbf{u})}\) with an appropriate smooth function \( \mathbf{W}: \mathbb{R}^d \times \mathcal{U} \subset \mathbb{R}^{n_u} \to \mathbb{R}^n\) . Along the unperturbed SSM \( \mathcal{A}_0\) , we then have

The dynamics on the perturbed aSSM \( \mathcal{A}_{\epsilon}\) can be written as

At leading-order, the aSSM is anchored along the slow trajectory \( \mathbf{x}^*(\epsilon t) \approx \mathbf{S}(\mathbf{u}(\epsilon t))\) which lies on the critical manifold \( \mathcal{L}_0\) . All trajectories of the aSSM-reduced system (13) will approach \( \mathcal{L}_0\) in forward time for any user-defined slow input \( \mathbf{u}(\epsilon t)\) with a small enough \( \epsilon > 0\) (see Fig. 1b). Note that we have chosen the reduced coordinate such that the critical manifold \( \mathcal{L}_0\) satisfies \( \mathbf{r}=0\) .

In dynamic trajectory tracking problems, we are interested in determining the control inputs required to move the robot along a prescribed track. In contrast, in Section 2.2, we described a reduction methodology for the inverse problem in which a reduced model predicts a track under given control inputs. To solve the original trajectory tracking problem via this reduction, we can substitute the aSSM-reduced model for the full nonlinear dynamical system in the optimization problem (1). In that case, however, we need to employ an additional fast control component that seeks to eliminate the error of the slow, aSSM-based controller.

Specifically, we express the general control input as \( \mathbf{u}(t) = \mathbf{u}^s(\epsilon t) + \mathbf{u}^d(t)\) , where \( \mathbf{u}^s(\epsilon t)\) is a slow control input and \( \mathbf{u}^d(t)\) is a fast control input deviation. We assume \( \mathbf{u}^d(t)\) to be uniformly bounded, i.e., let \( |\mathbf{u}^d(t)| \leq \delta\) for some \( \delta>0\) and for all \( t \geq 0\) . We extend the non-autonomous system (2) to an autonomous system by letting

When \( \mathbf{u}^d(t)\equiv 0\) , system (18) is equivalent to system (4), for which we have already obtained the aSSM-reduced model (13). We will next show how this aSSM-reduced model carries over to system (18) for \( \mathbf{u}^d(t)\not \equiv 0\) small enough.

The aSSM, \( \mathcal{A}_{\epsilon}\) , of system (4) is constructed by the results of [1] as a small perturbation of the normally hyperbolic invariant manifold \( \mathcal{A}_0\) of system (6). These results do not directly imply that \( \mathcal{A}_{\epsilon}\) is also normally hyperbolic, but can be strengthened to do so under further technical assumptions (see Remark 3.319 of [1]). To avoid introducing technicalities beyond the scope of this paper, we simply assume that \( \mathcal{A}_{\epsilon}\) is normally hyperbolic. Then, by Theorem 3.1 of [1], \( \mathcal{A}_{\epsilon}\) perturbs into a nearby, diffeomorphic invariant manifold \( \tilde{\mathcal{A}}_{\epsilon}\) in system (18) for small enough \( \delta >0\) uniform bound in \( |\mathbf{u}_d(t)|\) . We depict the perturbed aSSM geometry of \( \tilde{\mathcal{A}}_{\epsilon}\) in Fig. 1c.

The perturbed aSSM-reduced model on \( \tilde{\mathcal{A}}_{\epsilon}\) at leading order in \( |\mathbf{u}^d|\) and \( \epsilon\) is given by

where \( \mathbf{B}(\mathbf{u}^s(\epsilon t)) = \mathbf{V}^\mathrm{T}(\mathbf{u}^s(\epsilon t)) \mathbf{D}_{\mathbf{u}}\mathbf{g}(\mathbf{S}(\mathbf{u}^s(\epsilon t)),\mathbf{u}^s(\epsilon t))\) can be interpreted as the effective reduced linear control matrix, that varies along the slow control inputs. As formula (20) shows, the leading order correction to the aSSM parametrization can be ignored when the control deviations are small. Equation (20) provides us with an aSSM-reduced model that replaces the full system in the finite-time horizon optimization problem (2.1).

We now combine our aSSM-reduced models with the finite-time horizon optimization problem defined in eq.(1). Recall, that the finite horizon optimal control problem takes as inputs the initial condition, the planning horizon, the target trajectory, workspace constraints, and actuation space constraints. Since, the planning horizon in practice is fairly short, the aSSM in this interval can be approximated to be the SSM attached to the initial condition of the optimization problem.

Using notation from section 2.1 and also the chart map of the critical manifold \( \mathcal{L}_0\) from Remark 1, we define the aSSM-reduced MPC scheme for the planning horizon \( (t_j,t_{j+1}]\) as

Here \( \mathbf{u}^s_0\) is the static control input needed to generate the static configuration \( \mathbf{x}(t_j)\) . This is dependent on observations made in the workspace after successful application of the optimal control inputs in the previous horizon \( (t_{j-1},t_{j}]\) . The equations of motion constraint for the above optimization problem can be derived by setting \( \mathbf{u}^s(\epsilon t) = \mathbf{u}^s_0\) and \( \mathbf{u}^d(t) = \mathbf{u}(t) - \mathbf{u}^s_0\) in the perturbed aSSM-reduced model eq.(20). By setting \( \mathbf{u}^s(\epsilon t) = \mathbf{0}\) and \( \mathbf{u}^d(t) = \mathbf{u}(t)\) in eq.(26) we arrive at the zeroth-order aSSM-reduced MPC scheme (see Remark 5). A first-order aSSM-reduced MPC scheme is obtained by assuming that the chart, parametrization and reduced dynamics are the same for all possible static control inputs \( \mathbf{u}^s_0\) in the aSSM-reduced MPC scheme (see Remark 6).

The optimization problem (26) can be solved, for instance, using the GuSTO method of [41]. GuSTO approximates the non-convex optimization problem into sequential convex optimization problems. The method has already been used in the context of solving SSM based MPC schemes in [42].

We outline our methodology and the control approach in Fig.5. The following sections will elaborate in detail the learning methodology for aSSMs.

Pan/Zoom

SSMLearn takes as input trajectory data collected from unforced experiments of the system. The data comprises decaying trajectories starting from a few initial conditions that ultimately converge to the equilibrium of the uncontrolled system. From a spectrogram analysis of this data and observed symmetries of the underlying physical system, we estimate the minimal dimension \( d\) of the dominant slow SSM to which these trajectories converge. With probability one, trajectories will not lie on the slow SSM, but truncating them to their final sections will provide an accurate approximation to the dominant slow SSM. We can use the Takens embedding theorem for delay-embedded observables to deduce the minimal embedding dimension \( p=2d+1\) for the SSM to embed smoothly in the \( p-\) dimensional observable space. The SSMLearn algorithm chooses the observable space’s coordinate frame such that resting position of the robot is always at the origin.

The truncated decaying data is arranged into a snapshot matrix \( \mathbf{Y} \in \mathbb{R}^{p\times N}\) . Furthermore, we define \( N_{\mathrm{traj}}\) as the number of collected decaying trajectories, if we assume each trajectory has length \( n_{i}\) , then \( N = \sum_{i=1}^{N_{\mathrm{traj}}} n_i\) . Performing singular value decomposition (SVD) on \( \mathbf{Y}\) and collecting the \( d\) dominant singular vectors gives an approximation of the tangent space \( \mathbf{V}(\mathbf{0})\in \mathbb{R}^{p \times d}\) of the zeroth-order aSSM. The zeroth-order aSSM parametrization and the autonomous part of its reduced dynamics are then obtained by polynomial regression. Specifically, we fit the expressions

via least-square minimization to trajectory data to determine the parameter matrices \( \mathbf{W}^\star_{\mathbf{0}}\) and \( \mathbf{R}^\star_{\mathbf{0}}\) . Here \( \mathbf{r}^{1:n}\) is a column vector comprising all the monomials built out of the coordinate components of \( \mathbf{A}\) from order \( 1\) to \( n\) . Therefore, we have \( \mathbf{r}^{1:n} \in \mathbb{R}^{1 \times m_{n}}\) with \( m_{n} = \binom{d+n}{d} -1\) . \( \mathbf{W}^\star_{\mathbf{0}}\) denotes the optimal parametrization coefficient matrix of size \( \mathbb{R}^{p \times m_{n_w}}\) and \( \mathbf{R}^\star_{\mathbf{0}}\) is the optimal reduced dynamics coefficient matrix of size \( \mathbb{R}^{d \times m_{n_r}}\) . The orders \( n_r\) and \( n_w\) are chosen in a way to avoid overfitting to the training data. We use the normalized mean trajectory error (NMTE) on test predictions \( \mathbf{y}_{\mathrm{pred}}(t)\) to optimize these dimensions, where

For more detail, see [44]. To learn the effect of the control deviation term on the SSM-reduced dynamics (22), we generate controlled training data by subjecting the robot to random time-varying control inputs \( \mathbf{u}(t)\) . We arrange the controlled data and the random inputs as snapshot matrices \( \mathbf{Y}_{u} \in \mathbb{R}^{p \times N_{t}}\) and \( \mathbf{U} \in \mathbb{R}^{n_{u} \times N_{t}}\) , where \( N_{t}\) is the length of the controlled trajectory. The effective reduced linear control matrix \( \mathbf{B}(\mathbf{0})\) is then selected as the solution to the minimization problem,

Using this data-driven, zeroth-order aSSM-reduction methodology as a template, we will next discuss methods to learn a \( d\) -dimensional aSSM-reduced model’s chart map, parametrization and reduced dynamics.

3.2.1 Pointwise static SSM models

3.2.2 Control calibration for static SSMs

3.2.3 Interpolating an aSSM-reduced model from static SSM models

3.2.4 Systematic parameter selection for aSSM-reduced models

We apply the method outlined in Sections 2.2-2.3 to construct a perturbed aSSM-reduced model (20) for open-loop control of a double pendulum. We consider the observable space \( \mathbf{y} \in \mathbb{R}^{p}\) to be the full phase space variable \( \mathbf{x} = {\theta_1,\theta_2,\dot{\theta}_1,\dot{\theta}_2}\) , which renders \( p =4\) . We use the governing equations of the double pendulum derived by [45], who model the double pendulum as a set of two connected rods. We set the lengths and weights of the rods equal and to \( 1 \text{ [m]}\) and \( 1 \text{ [kg]}\) , respectively. The system is subject to a gravitational acceleration of \( g = 9.8  [\text{m/s}^2 ]\) . We assume linear Rayleigh damping for each rod, with the top rod damped twice as much as the lower one.

The dimension of the control input \( \mathbf{u}\) is \( n_u = 2\) , representing the torques acting on each angular degree of freedom of the system. We choose the grid coordinates for interpolation to be the angular coordinates of the rods \( \mathbf{q} = {\theta_1,\theta_2}\) . We also perform a linear stability analysis across static torque configurations and find a square region \( \theta_1,\theta_2 \in [-\frac{\pi}{9},\frac{\pi}{9}]\) that corresponds to stable static states. We also observe that the assumption in eq.(9) is satisfied in this region.

We generate two, \( 100\) -second long decaying training trajectories for \( N_u=36\) static torque configurations using random initial conditions. The static states live on a square grid \( \mathbf{Q}_s \in [-\frac{\pi}{9},\frac{\pi}{9}]^2.\) We observe that after a suitable truncation of the dataset, a \( d=2\) dimensional SSM-reduced model captures the decaying dynamics accurately for all static torques. We identify static states with \( \theta_1= -\theta_2\) tend to have more nonlinear SSMs compared to others. Setting \( n_w = 3\) and \( n_r = 3\) provides us with good static SSM approximations that have average \( \mathrm{NMTE} \approx 4\% \) across all test data. Since our interpolation grid is square and regular, we employ a spline-based MATLAB interpolation for the aSSM-reduced model’s coefficients. For the control deviation term \( \mathbf{u}_d(t)\) in eq.(20), we can derive the exact form of the linear control matrix as \( \mathbf{B}(\mathbf{q}) = \begin{pmatrix} \mathbf{0} \\ \mathbf{V}^\mathrm{T}(\mathbf{q})\end{pmatrix} \in \mathbb{R}^{4\times 4}\) , because we use the phase space as observable space.

We now compare our learned perturbed \( 2\mathrm{D}\) aSSM-reduced model with the full model simulation for slow control torques that generate a random slow trajectory in \( \mathbf{Q}_s\) , combined with a chaotic control deviation. The random slow trajectory \( \mathbf{S}(\mathbf{u}^s(\epsilon t))\) is generated by a gradient noise scheme called Perlin noise (see [46]) and is rescaled to fit into the grid window \( \mathbf{Q}_s\) . The chaotic control deviation is \( \mathbf{u}^d(t) = \delta \tilde{x}(t)\) , where \( \tilde{x}(t)\) is a unit rescaled scalar trajectory of one component from the response of a \( 3\mathrm{D}\) chaotic Lorenz model (see [47]) launched from a random initial condition.

Pan/Zoom

In Fig. 6a, we compare the different orders of aSSM approximation for an initial condition launched on the aSSM for \( \delta =0\) and \( \epsilon = 0.005\) . The aSSM and its first-order approximation provide high-accuracy predictions. The zeroth-order approximation does well initially but loses accuracy in the long term. Figure 6b depicts the case when \( \delta=1.2\) , showing that the perturbed aSSM-reduced model and its first-order approximation accurately match the full model. Selecting \( \delta=3.2\) yields a similar result, thus proving the robustness of aSSM-reduced models (see Fig. 6c). We also plot in Figs. 6d-e snapshots of the \( 2\mathrm{D}\) aSSM slices and dynamics in the full phase space. These plots highlight the geometric significance of the aSSM, showing that the true solution lies close to the aSSM slice and is accurately captured by the aSSM-reduced dynamics. The aSSM slice is anchored to the slow target and translates along it changing shape, whereas the prediction of the zeroth-order aSSM approximation remains close to a fixed aSSM slice for all times.

We collect \( 3\) -second long samples of \( 40\) decaying trajectories for each of \( N_u=100\) statically forced equilibrium configurations of the trunk. We split the \( 40\) decaying trajectories into \( 15\) for training and the rest for testing. This gives a total of \( 1500\) training trajectories for the whole workspace of the robot. Each of the \( 100\) static control input vectors \( \mathbf{u}^s\) was generated by random sampling of \( 100\) distinct combinations of \( 8\) values from the set \( {0,150,350}\) . \( 75\) of these static steady states densely cover a spherical cap having radius \( 20 \text{ [mm]}\) and height \( 5 \text{ [mm]}\) in the workspace (see Fig. 8a). We also ensure the combination corresponding to the zero control input vector is part of this collection.

We concentrate on the collected trunk decay data about the unforced configuration. We first split each set of \( 40\) collected trajectories into \( 15\) for training and the rest for testing. Next, we truncate the training data to half a second. We follow this up with the construction of a \( 3p\) - dimensional delay embedded observable space \( y = (\mathbf{z}(t),\mathbf{z}(t+\tau), \dots ,\mathbf{z}(t+(p-1)\tau))^\top \in \mathbb{R}^{3p}\) , where \( \tau = 0.01 \text{ [s]}\) . We observe that for \( p=4\) , a \( d=5\) dimensional SSM-reduced model for \( n_w = 4\) and \( n_r = 2\) accurately captures the decay in a \( 12\) -dimensional delay embedded observable space, with an average \( \mathrm{NMTE} \approx 3 \%\) on test data (see Fig. 8b).

Pan/Zoom

Since the \( z_{ee}\) coordinate is directly affected by gravity, the decay data in the coordinate is over-damped. This is reflected by a real eigenvalue of the linear part of our 5D SSM-reduced model. Specifically, the spectrum of the linear part of the \( 5\mathrm{D}\) zeroth-order aSSM approximation is

The eigenvalues \( \lambda_1\) and \( \lambda_3\) are not in exact resonance, but the gap between them is small compared to \( \lambda_5\) . This near-resonance arises from the \( x_{ee}-y_{ee}\) symmetry of the trunk’s equilibrium configuration. The resonance is not exact because the finite-element model of the robot is not c onstrained to respect this symmetry.

If we truncate the training trajectories at \( 1\) second, we find that a \( 4\) D zeroth-order aSSM approximation already captures the local dynamics accurately. The spectrum of the reduced dynamics in that case is

The truncation allows us to get a refined spectrum that is near resonant up to the first decimal but comes with the cost of losing the overdamped signatures of the trunk.

Using the same parameters for the learned 4D and 5D zeroth-order aSSM approximations, we construct the static SSM models for 4D and 5D aSSM-reduced models. We will compare these models in open-loop tests and choose the best performing one for closed-loop tasks.

We further collect controlled trajectory data following the procedure outlined in Section 3.2.2. The random control deviation \( \mathbf{u}^d(t)\) , with a maximum norm of \( \delta = 50\) , is generated using the Latin hypercube sampling technique (see [49]). We generate one 200-second long random control input deviation. To this we add the random 100 static inputs from earlier and pass the resulting control inputs to the robot to generate controlled responses. We use these controlled responses and their corresponding inputs to learn local linear control matrices \( \mathbf{B}_s\) by solving eq.(34). We repeat the same procedure but we keep the autonomous reduced dynamics fixed to the zeroth-order aSSM approximation to arrive at \( \mathbf{B}_s\) for the first-order aSSM-reduced model.

For learning \( \mathbf{B}_{\mathbf{0}}\) (see eq.(31)) for the zeroth order aSSM-reduced approximation, we generate a random control deviation \( \mathbf{u}^d(t)\) with \( \delta = 350\) . We train this model for larger magnitudes of inputs to better anticipate large control deviations that will arise in our closed loop control tasks.

To test our aSSM-reduced models on unseen general control inputs, we generate new random control input sequences \( \mathbf{u}(t)\) with a maximum norm of \( \delta = 350\) . The trajectory generated by this input traces out a spherical cap in the workspace. We further generate faster versions of the same random input \( \mathbf{u}(t)\) that result in faster mean speeds achieved by the trunk robot in the workspace. We run controlled simulations with \( \mathbf{u}(t)\) and its faster versions on the trunk robot.

On each controlled response and their corresponding inputs, we perform a detailed validation exercise given by :

• Randomly select \( 50 \%\) points from the control response trajectory.
• Use these points as initial conditions to construct an aSSM-reduced model appearing in eq.(26). The construction of the aSSM-reduced model is based on interpolation or regression methods discussed in Section 3.2.3.
• Simulate this model for 0.05 seconds using the corresponding general control input values that occur within the time window of the actual response, starting from the same initial condition.
• Calculate the NMTE error metric between the true response of the trunk and the prediction from the aSSM-reduced model.
• Plot the initial conditions and shade the plot with a color gradient, where low NMTE values are represented by green and high NMTE values by red.

Further, we perform the above validation on the following aSSM-reduced interpolation models :

• \( 5\) D aSSM-reduced model with modified inverse distance weighting (MIDW) interpolation, with \( p=2\) in eq.(35), constructed on a scattered grid in workspace coordinates \( \mathbf{q} =\mathbf{z} \in \mathbb{R}^3\) .
• \( 5\) D aSSM-reduced model with quadratic polynomial regression (QPR), setting \( n_q =2\) in eq.(36), constructed on a scattered grid in the first two reduced coordinates \( \mathbf{q} ={(\mathbf{V}^{\top}_0\mathbf{y})_{1},(\mathbf{V}^{\top}_0\mathbf{y})_{2}} \in \mathbb{R}^2\) .

We compare these results with the \( 4\) D aSSM-reduced model counterparts and their zeroth-order approximations as well.

Pan/Zoom

Figure 9 shows NMTE color plots of the performance of the \( 4\) D models. Notice that these models perform poorly for initial conditions having a non-zero \( z_{ee}\) value. This is expected as these models do not have predictive power that models the overdamped behavior of the \( z_{ee}\) direction. In Fig. 10, we have the results for the \( 5\) D aSSM-reduced models. In comparison to the \( 4\) D case, we see nearly complete green patches signaling low NMTE values for nearly all the initial conditions. This informs us that the our \( 5\) D aSSM-reduced models have significantly greater generalizability compared to the \( 4\) D aSSM-reduced models. At the same time, the plots signal that \( 4\) D aSSM-reduced model is an optimal choice for closed-loop control tasks restricted to the \( x_{ee}-y_{ee}\) plane.

Pan/Zoom

When comparing the MIDW interpolation and QPR methods, we observe that MIDW offers slightly better performance than QPR for the \( 5\mathrm{D}\) aSSM-reduced model. However, this comes with longer mean run time for MIDW due to complexity in the interpolation procedure. The QPR method has lower mean run time which aids in achieving closed-loop performance in real-time, and it also comes with good prediction accuracy on unseen inputs. Based on these findings, we choose the \( 5\) D aSSM-reduced model approximated using QPR in our MPC schemes for closed loop control.

The zeroth-order aSSM-reduced model’s performance also meets our expectation. Indeed, a small region close to the origin has green patches in Fig.10 but as we move further away from the origin the errors start to grow.

For the above results, the average instantaneous speed of the trunk’s tip during the full control response is approximately \( \approx 125 \text{ [mm/s]}\) . This gives the non-dimensional slowness measure \( r_s \approx 1.7\) . We show similar results when the average instantaneous speed is about \( 265 \text{ [mm/s]} \) with an \( r_s \approx 3.7\) value for the \( 5\) D aSSM-reduced models in Appendix B.

Our first target \( \boldsymbol{\Gamma}(t)\) is a figure-8 trajectory in the \( x_{ee}-y_{ee}\) plane. This figure-\( 8\) target has a 10-second period and a maximum amplitude of \( 40 \text{ [mm]}\) . The average instantaneous speed achieved along the target is \( 37.74 [\text{mm/s}^2]\) with slowness measure \( r_s \approx 0.5\) . The target trajectory starts at the origin and traverses the figure-\( 8\) in the clockwise direction. Figure 11 shows results from four different models used in a MPC control scheme: \( 5\) D aSSM-reduced model, \( 4\) D aSSM-reduced model, \( 5\) D zeroth-order aSSM-reduced model and \( 5\) D first-order aSSM-reduced model. We find that the \( 4\) D aSSM-reduced model has the lowest ISE. Using this as a baseline, we evaluate relative ISE for the rest of the models (see Fig. 11b). These results are consistent with the open-loop analysis presented in Section 4.3 in terms of the performance of the various models.

Despite the relative ISE of \( 206 \%\) , the \( 5\) D aSSM-reduced model also achieves similar control performance along the target as the \( 4\) D case. We also observe that the \( 5\) D first-order aSSM-reduced model does better than the \( 5\) D zeroth-order aSSM-reduced model. The \( 5\) D zeroth-order aSSM-reduced model traces the target near the origin with good accuracy, but as it moves radially outward, it starts to overestimate the target.

Pan/Zoom

On the other hand, the first-order aSSM-reduced model tracks the majority of the figure-8 precisely but begins to jitter on the bottom curve on the right lobe causing the relative ISE to shoot up to \( 2,801 \%\) . We attribute this phenomenon to model error arising from the control calibration performed in a region in the workspace where the first-order aSSM approximation is no longer valid.

We conclude aSSM-reduced models are generalizable, as none of these models were trained on inputs that resembled the figure-\( 8\) . We also find that our aSSM-reduced models to extrapolate well in regions of the \( x_{ee}-y_{ee}\) plane, \( \sqrt{x_{ee}^2 +y_{ee}^2}>20 \text{ [mm]}\) , where they have no prior information of the static steady states and SSM computations. In Appendix E, we present comparisons of the aSSM-reduced models with linear methods. The results from Appendix E indicate that the trunk is only mildly nonlinear in \( x_{ee}-y_{ee}\) . Indeed, both linear methods are able to trace out the figure-\( 8\) but lack precision compared to the aSSM-reduced models.

Pan/Zoom

We also investigate the trade-off between model accuracy and finite-horizon optimization solve times. Figure 12 plots the Pareto front for the figure-\( 8\) trajectory. The Pareto front here, qualitatively depicts the trade-off between control accuracy and control solve times. The models appearing on the Pareto front are said to be Pareto optimal, as their performance in terms of accuracy and solve time is not strictly dominated by any other model. We find the \( 4\) D aSSM-reduced model and the \( 5\) D zeroth-order aSSM-reduced model to lie on the Pareto front. The \( 4\) D aSSM-reduced model provides us precision control with average solve time of about \( 31 \text{ [ms]}\) whereas the \( 5\) D zeroth-order aSSM-reduced model in comparison has a \( 10 \text{ [ms]}\) lower solve time but comes at a loss of control accuracy as it has the largest ISE. Depending upon the use case for the robot, the practitioner can make a choice between the two Pareto efficient models for planar trajectory tracking.

We now consider a \( 3\) D target track in the workspace, randomly generated using a Perlin noise scheme (see [46]). The track starts at the origin and varies smoothly in all workspace directions. The target \( \boldsymbol{\Gamma} = {\alpha(t),\beta(t), \sqrt{R^2 - s^2(\alpha(t)^2 - \beta(t)^2)} }\) , where \( \alpha(t)\) and \( \beta(t)\) are unit normalized scalar signals, generated by such a scheme. We set the maximal radial span of the target at \( s=40 \text{ [mm]}\) and the maximal height at \( h = 20 \text{ [mm]} \) , which gives \( R=\frac{h^2+s^2}{2h}\) .

To test our aSSM-reduced MPC schemes for speed variations along the given track, we consider a mean instantaneous speed of \( 50.16 [\text{mm/s}^2]\) with slowness measure \( r_s \approx 0.7\) and faster variation at a mean instantaneous speed of \( 100.32 [\text{mm/s}^2]\) with \( r_s \approx 1.4\) . In Fig. 13a, we show closed-loop performance for the slower varying case using the aSSM-reduced models as well as linear models. We find the aSSM-reduced model and its approximations to perform significantly better than the linear methods. TPWL method performs visibly poorly and is unable to reach the maximal height. The Koopman static pregain method is able to roughly trace out the random target’s behavior but fails to follow it accurately.

Pan/Zoom

Pan/Zoom

Figure 14a confirms our visual assessment quantitatively with the relative ISE bar plots indicating the same trend. Results for the faster version of the same track are presented in Fig. 13b. Again, a visual assessment suggests that the aSSM-reduced model has the best performance, confirming the robustness of our aSSM-reduced model. We again verify this assessment quantitatively from the ISE bar plots for the various models in Fig. 14b. As seen from the plot, Koopman static pregain does significantly worse here than in the slow case. This is because the method assumes the inputs to be quasi-statically varying and the mapping between the static inputs and steady states to be linear, neither of which holds for the trunk. Generally, the key reason for the failure of the linear methods is their non-robustness to unseen control inputs.

We compare more closely the aSSM-reduced models and their approximations by plotting the performance in the \( z_{ee}\) direction in Fig. 15. The results confirm the local domain of validity of the \( 5\) D zeroth-order aSSM-reduced model showing degraded performance along the height of the workspace for both speeds of the target track. The first-order aSSM-reduced model shows a slight improvement in performance but fails to reach the heights of the target, even though it tracks the variations smoothly. We again note that \( z_{ee}\) values explored by the target are much further away from the training regime of the static SSM models.

In Fig. 16, we plot the Pareto fronts for the slow and fast random tracks. Similarly to the figure-8 results, we find the zeroth-order aSSM-reduced and the aSSM-reduced model to be Pareto efficient. We also notice the average solve times of the aSSM-reduced model decreases with increasing target speed but the ISE increases. Note that the faster target generates a lower ISE under the zeroth-order aSSM approximation. This is because this approximation is robust to uniformly bounded variations in the control inputs. The first-order aSSM approximation is suboptimal because for a slight improvement in ISE, the average solve time is still larger compared to the aSSM-reduced model.

Pan/Zoom

Pan/Zoom

In the next example, we increase the difficulty of the control task with a \( 3\) D Pacman track that is mostly outside the training regime of all the models and also has non-smooth variations. The radius of the Pacman track is \( 30 \text{ [mm]}\) and its height is \( 25 \text{ [mm]}\) . The mean instantaneous speed along the track is \( 22.9 \text{ [mm/s]}\) , which gives the slowness measure \( r_s \approx 0.3\) .

Pan/Zoom

Figure 17a shows the closed-loop performance of the aSSM-reduced models. A visual assessment shows that the aSSM-reduced model performs the best, as expected. The other two approximations struggle to reach the maximum height, although the first-order approximation performs better than the zeroth-order. The first-order aSSM-reduced model is at least able to track the steep slope, but ultimately fails to stay on the circular track. The first-order aSSM approximation results illustrate that even if one knows the existence and approximate location of an aSSM, one generally cannot approximate it as a fiber bundle of constant linear fibers on larger domains because aSSMs tend to be fiber bundles of non-constant nonlinear fibers.

In Fig. 17b, the relative ISE plots confirm our visual assessment. For completeness, in Fig. 18, we also plot the \( z_{ee}\) tracking direction, which yields similar conclusions. The Pareto front is again similar as in the previous examples (see Fig. 19). The aSSM-reduced model has the lowest ISE but has a slightly increased mean solve time of \( 0.32 \text{ [s]}\) compared to the zeroth-order aSSM at \( 0.24 \text{ [s]}\) .

Pan/Zoom

Pan/Zoom

As a final example, we consider a similarly shaped Pacman track with the same radius but at a height of \( 15 \text{ [mm]}\) and with five randomly sized and placed spherical constraints in the workspace. The average speed clocked along this track is \( 21.8 [\text{mm/s}^2]\) with slowness measure \( r_s \approx 0.3\) . The goal is to control the trunk as closely as possible along the track but at the same time avoid all the constraints.

Pan/Zoom

In Fig. 20a, we plot the closed loop predictions for the aSSM-reduced model and its approximations. We also plot the spherical constraints as red transparent spheres in the workspace. We introduce two error metrics that will aid in our quantification. The first is violation ratio, which is the number of times the robot entered the constraint divided by the length of the robot’s trajectory. The second is maximum violation, which is the maximum distance by which the robot violates the constraints. In Figs. 20c-d, we plot these metrics for the aSSM-reduced model, the first-order aSSM-reduced model and the zeroth-order aSSM-reduced model. The radii of the spherical constraints range from \( 1 \text{ [mm]}\) to \( 5 \text{ [mm]}\) . By maximal violation measure, all the aSSM-reduced models perform consistently; as the maximum violation distance is \( 0.32 \text{ [mm]}\) across models. We remark that these are indeed hard constraints, so in principle, when they are violated, the program algorithm should exit. However, to test the models across the target horizon, we re-initialize the failed models at an initial condition just outside the constraint, closest to the point of failure.

Overall, the aSSM-reduced model has the smallest violation ratio, the smallest maximum violation distance of \( 0.01 \text{ [mm]}\) and also the lowest ISE. This effectiveness of the aSSM-reduced model for constrainted workspace configurations illustrates that the mathematical assumptions going into the theory of aSSM-reduction are realistic, at the very least for soft trunk robots.

Our focus is to control the arm’s free end (end effector), hence our workspace is defined by the end effector coordinates \( \mathbf{z} = {x_{ee}, y_{ee}, z_{ee}} \in \mathbb{R}^3\) . We collect \( 10\) -second long samples of 4 decaying trajectories for each \( N_u = 1000\) randomly sampled static arm configurations. We ensure practical limitations on observability by only assuming access to a \( 3p\) -dimensional observable space \( y = (\mathbf{z}(t),\mathbf{z}(t+\tau), \dots ,\mathbf{z}(t+(p-1)\tau))^\top \in \mathbb{R}^{3p}\) , where \( \tau = 0.0168 \text{ [s]}\) .

The undeformed arm is symmetric along the \( x_{ee}-y_{ee}\) and also freely stretches along the \( z_{ee}\) direction, due to these reasons we expect three dominant slow modes in the decay profiles of the arm. By Takens’ delay embedding theorem, we require \( 3p \geq 15\) to guarantee capturing decay along these three dominant slow modes. Indeed for \( p=5\) , we observe \( d=6\) dimensional static SSM-reduced models with \( n_r = 2\) and \( n_w = 2\) to accurately capture the decay with an average \( \mathrm{NMTE} = 5 \%\) on test data for the short arm and an average \( \mathrm{NMTE} = 4.5 \%\) on test data for the long arm .

We further collect controlled data by adding a 10-second long randomly generated control deviation input, with a maximum torque magnitude of 1.5, to a step input sequence formed by 100 static inputs, each held for 10 seconds (see Section 3.2.2 for details). We use this to learn the effective control calibration set. Specifically, we find a \( 6\) D aSSM-reduced model obtained via quadratic polynomial regression (QPR) on the static SSM coefficient set to provide consistent performance on unseen test data. For a 30-second long random input with slowness measure \( r_s = 0.75\) , our \( 6\) D aSSM-reduced model’s prediction has a \( 9 \%\) NMTE error for the short arm and a \( 11 \%\) NMTE error for the long arm. See Appendix F for prediction results in the robot’s workspace.

We design a spherical target track \( \boldsymbol{\Gamma} = {r(t)\cos(2\pi t/T),r(t)\sin(2\pi t/T), \sqrt{1 - r^2(t)} }\) , where \( r(t) = r_0 (1+\cos(6\pi t/T + \pi))\) . This track resembles a three-lobed curve, known as the Trifolium curve. We set the Trifolium curve’s parameters \( r_0=0.45 \text{ [m]}\) and \( T=10 \text{ [s]}\) . The mean instantaneous speed along the track is \( 56.8 \text{ [cm/s]}\) and the slowness measure is \( r_s =0.5\) .

We perform closed-loop control using the reduced MPC scheme eq.(26) for the 6D aSSM-reduced model and its first-order approximation. We also add a comparison with the Koopman static pregain method, following similar guidelines in Appendix C.2 using identical training data collected for the static SSM models. The planning horizon for the MPC scheme is set to \( 0.1 \text{ [s]}\) and the workspace cost matrix \( \mathbf{Q}_z = 10 \cdot\mathbf{1}_{3 \times 3}\) for all models. We set \( \mathbf{R}_u = 0.8\cdot\mathbf{1}_{6 \times 6} \) for the 6D aSSM-reduced MPC scheme and \( \mathbf{R}_u = 1.5\cdot\mathbf{1}_{6 \times 6} \) for its first-order approximation and the Koopman static pregain method. For this example, we observed that different \( \mathbf{R}_u\) values led to improved performance for the first-order aSSM-reduced approximation and the Koopman static pregain method; hence, the \( \mathbf{R}_u\) settings were adjusted to ensure a fair comparison.

In Fig.21a-c, we plot the closed loop predictions in the robot’s workspace for all the model-based MPC schemes. It is impossible for the soft arm to exactly track the spherical Trifolium curve as the soft arm is underactuated (see Appendix F). Hence, all plots show a gray-shaded cylindrical tube around the target, indicating the allowed operational tolerance of the soft arm’s end effector during closed-loop tracking. We set the radius of the operational tolerance tube to be \( 2\%\) of the target’s length, which amounts to a tube of radius \( 15 \text{ [cm]}\) . We clearly observe that only the 6D aSSM-reduced model falls within the operational tolerance tube. Both the 6D first-order aSSM-reduced model and the Koopman static pregain method produce closed-loop inputs that cause the arm to deviate from the target trajectory as it approaches the maximal points.

Pan/Zoom

Pan/Zoom

Pan/Zoom

We support these observations in Fig.25a-b, by plotting bar plots for the relative ISE and the Pareto plots comparing ISE with solve times. Across all metrics, the 6D aSSM-reduced model offers the lowest ISE \( \approx 3 [\text{m}^2 s]\) and the lowest average solve time \( \approx 1 \text{ [ms]}\) . The first-order aSSM-reduced approximation has comparable solve times but overall accumulates errors, leading to larger ISE due to the model’s inaccuracy in the extended target regions. The linear Koopman static pregain tracks worse in the extended target regions because the linear approximation is invalid there. Specifically, the linear Koopman model’s limitations become evident when it extrapolates to large control inputs, leading to increasing instability in the soft arm’s response. The controller counteracts this instability by clipping the predicted control inputs (see Fig.21c). In contrast, our results show that nonlinear aSSM-reduced models achieve robust extrapolation accuracy in these regions (see Fig.21a-b), relying solely on their inherent nonlinearity to predict optimal control inputs without requiring corrective control measures. In Fig.28a-c, we also present snapshots of the elastic arm tracking for each model’s prediction and display the accumulated normalized ISE The ISE is normalized by dividing it by the target’s total time and the square of the undeformed arm’s length and is denoted as \( \overline{\mathrm{ISE}}\) ..

Pan/Zoom

Pan/Zoom

Pan/Zoom

Pan/Zoom

Pan/Zoom

We have not included comparisons with the zeroth-order aSSM approximation, as we know from theory that the method is inapplicable to track large target tracks, which we have systematically demonstrated for the soft trunk robot example in Section 5. We also note that the TPWL method of [50] cannot be implemented in the elastic arm setup, as it relies on a finite-element model of the arm, which is not readily available.

Our closed-loop MPC scheme for the elastic arm uses optimized evaluation techniques for the aSSM-reduced model, yielding a 6D aSSM-reduced model with significantly lower average solve times of 1 [ms] for a 0.06 [s] MPC horizon. These solve times are 10 times lower than those reported in [37] for hardware MPC experiments with simpler 6D SSM-reduced models. With these upgrades, higher-dimensional aSSM-reduced models offer better accuracy with faster solve times. We note that shorter solve times alone do not necessarily translate into improved performance in high-speed tasks. Indeed, the robot’s sampling time is typically hardware-constrained irrespective of the specific task. Therefore, what ultimately matters in a hardware implementation is the reduced model’s accuracy in predicting the dynamics.

In a hardware setting, the elastic arm’s end effector is tracked by a motion capture system, which approximates its position as a rigid body using markers embedded on the elastic arm. OptiTrack systems are commonly used in academic robotic hardware systems to record observations. They offer millimeter-accuracy when there are enough cameras. Their placement ensures successful calibration without occlusions and there are no external light sources that introduce noise. However, in most practical hardware applications, it is difficult to meet these requirements.

From the OptiTrack documentation [53], a poorly calibrated setup results in an error range of approximately 0.5-1 [cm]. To mimic this practical limitation in our simulations, we add noise sampled from a bounded Gaussian distribution with an upper bound of 1 [cm] to the feedback and use it as the initial condition for the MPC. To account for this additional operational uncertainty, we increase the operational tolerance radius by 1 centimeter.

Pan/Zoom

Pan/Zoom

Pan/Zoom

Pan/Zoom

We test the closed-loop performance of all three models for the Trifolium target track discussed in the previous section in the presence of feedback noise discussed above. In Fig. 32a-c, we plot the closed-loop predictions in the robot’s workspace. We observe that camera feedback noise causes the short elastic arm to jitter in the 6D aSSM-reduced model and in its first-order approximation. The jitter is not visible in the Koopman static pregain model, as the method is dominated by the linear mapping that acts directly on the target track, and the dynamic feedback contribution is minimal and does not significantly affect it. Overall, the 6D aSSM-reduced model remains within the operational tolerance tube when the other methods tend to exit it and develop larger ISE errors (see Fig. 32d).

We conclude this section with a closed-loop demonstration for the long elastic arm on a spiral-shaped target track with a mean instantaneous speed of \( 103.3 \text{ [cm/s]}\) and slowness measure \( r_s = 1.3\) . The track starts from the undeformed arm’s configuration and spirals downward to a height of \( 2.3 \text{ [m]}\) . The spiral track is depicted as a dotted curve in Fig. 37a. For this case, we perform closed-loop comparisons for the 6D aSSM-reduced model and the Koopman pregain static method. We set the planning horizon for the MPC scheme as \( 0.1 \text{ [s]}\) and the workspace cost matrix \( \mathbf{Q}_z = 100 \cdot\mathbf{1}_{3 \times 3}\) for both models. We identified that \( \mathbf{R}_u = 0.5\cdot\mathbf{1}_{8 \times 8} \) for the 6D aSSM-reduced MPC scheme and \( \mathbf{R}_u = 0.01\cdot\mathbf{1}_{8 \times 8} \) for the Koopman static pregain method produced reasonable optimal control inputs that led to the lowest errors during closed-loop. We omit the first-order approximation of the aSSM-reduced model here because the results showed performance gains similar to those already demonstrated for the short elastic arm.

Pan/Zoom

Pan/Zoom

Pan/Zoom

Pan/Zoom

Pan/Zoom

Figure 37a-b shows snapshots of the closed-loop results of the long elastic arm for the 6D aSSM-reduced model and the Koopman static pregain method. In both plots, we denote the operational tolerance tube (OTT) in gray, with a radius of 17.5 [cm] amounting to \( 2\%\) of the spiral target’s length. We observe that the aSSM-reduced model stays within the OTT and provides the lowest ISE (see Fig. 37c-e) compared to the Koopman static pregain method. The example yet again demonstrates the limitations of linear methods for accurately tracking the target in extended regions, where the trunk assumes nonlinear shape profiles to course-correct and achieve lower errors, which our nonlinear aSSM-reduced model can do. In Appendix F.3, we also briefly present the long elastic arm’s closed-loop results in the presence of experimental noise, as in the short elastic arm’s case, our aSSM-reduced model robustly tracks the target track.

he soft trunk robot geometry is densely meshed using smaller polyhedra or hexahedra 3D objects. We denote the generalized displacement of all nodes in the meshed trunk description as \( \mathbf{n} \in \mathbb{R}^{n_f}\) . The mesh material needs to be deformable and is governed by a nonlinear elastic continuum law that depends on the material’s Young’s modulus and Poisson’s ratio, determined from the linear elastic regime.

Along certain mesh regions, nodes make contact with cables actuators, the forces imparted on the mesh due to cable actuation are given by the function \( \mathbf{B}(\mathbf{n}) \mathbf{u}\) , where \( \mathbf{u} \in \mathbb{R}^{n_u}\) denotes the control input vector. Employing the finite element method with these mesh properties and cable actuation constraints, one can approximate the partial differential equations describing the soft trunk into a finite-dimensional second-order ordinary differential equation,

Here \( \mathbf{M}(\mathbf{n})\) is the mass matrix, \( \mathbf{F}(\mathbf{n},\dot{\mathbf{n}})\) is the nonlinear internal forcing, \( \mathbf{B}(\mathbf{n}) \mathbf{u}(t)\) is the forcing due to cable actuation, and \( \mathbf{F}_{ext}(t)\) is a generalized external forcing acting on the soft trunk. In our trunk setup, the mesh consists of \( n_f = 709\) nodes and \( n_u = 8\) cables, with the external force held constant at the Earth’s gravitational force. For exact details on the implementation of the nonlinear elastic internal force and the actuation force, see [48].

The soft arm is modeled as a Cosserat rod with undeformed length \( l\) , cross-sectional area \( A\) , mass per unit length \( \rho\) , second area moment of inertia \( \mathbf{I}\) , and bending and shearing stiffness matrices \( \mathbf{B}\) and \( \mathbf{S}\) .

We denote the centerline coordinate of the rod as \( \mathbf{x}(s,t) \in \mathbb{R}^3\) , the angular momentum of the rod \( \boldsymbol{\omega}(s,t) \in \mathbb{R}^3\) , and the generalized curvature of the rod \( \boldsymbol{\kappa}(s,t) \in \mathbb{R}^3\) , where \( s \in [0,l]\) is a material coordinate and \( t\) is time. We also define the rotation matrix \( \mathbf{T} \in \mathbb{R}^{3\times 3}\) that transforms vectors from the lab frame to the rod’s local frame. Then, the governing PDEs of the soft arm are given by

where \( \hat{\boldsymbol{\omega}}\) is the skew-symmetric matrix of \( \boldsymbol{\omega}\) , \( \mathbf{F}_{SS}\) represents shearing or stretching forces, \( \mathbf{F}_{B}\) represents bending or twisting forces, \( \mathbf{F}_{SC}\) is an internal coupling force, \( \mathbf{F}_{AD}\) represents advection and dilation forces, \( F_{ext}\) is the external force and \( C_{ext}\) is the external torque. To mimic realistic damping in these systems, a proportional damping term per unit length is introduced after spatial discretization of the soft arm’s PDE. The external torques act on user-defined control points along the rod’s length, which serve as the control actuation forces for the arm. For exact forms of the internal force terms and the damping force terms, refer to [52]. For the actuation settings related to our soft arm example, refer to [15].

\color{black}

We closely follow the TPWL method implemented in [50]. We outline the steps as follows:

• Step 1: Collect 5 random control input sequences with a maximum value \( \delta = 250\) using a Latin hypercube sampling technique. Subject the trunk robot to these input values.
• Step 2: Set up a snapshot matrix for the controlled responses of the trunk. Perform SVD on this dataset and collect singular vectors corresponding to the \( 28\) largest singular values and form the reduced basis.
• Step 3: Evaluate the Jacobian of the full model at the initial condition of the random control responses. Reduce the Jacobian using the learned basis from Step 2.
• Step 4: Run the SVD-reduced linear model subjected to general control inputs and compare the predictions with the full response.
• Step 5: Stop step 4 when the errors in the predictions exceed a user-defined threshold. Store the SVD-reduced linear model reduced by SVD and the initial condition obtained from step 3.
• Step 6: Repeat steps 3 and 4 at the violation point. Stop after scanning all the random control response trajectories.

We evaluate the error as a pointwise relative error for both velocities and positions and set the threshold to \( 10 \%\) . Here, our method differs from the original TPWL method, which only evaluates the errors in velocities for a threshold value of \( 200\) . Our settings resulted in \( 52\) piecewise \( 28\) -dimensional linear models.

We implement the Koopman static pregain method described in [27] for the bionic trunk robot. We choose the observable space to be a delay embedding of the end effector position \( \mathbf{y}={\mathbf{z}(t),\mathbf{z}(t+\tau)} \in \mathbb{R}^{6}\) , with \( \tau = 0.01\) . We generate controlled training data by passing a step input signal that combines all static input probing signals used to train the aSSM-reduced model. We arrange the trajectories and the control input as snapshot matrices \( \mathbf{Y} \in \mathbb{R}^{6 \times N_t}\) and \( \mathbf{U}\in \mathbb{R}^{n_u \times N_t}\) , respectively, where \( N_t\) is the total length of both trajectories and \( n_u\) is dimension of the robot’s actuation space. We fit a linear time invariant system to this dataset:

This is defined as the dynamic Koopman operator which is a linear non-autonomous system in the observable space. A random controlled trajectory will behave linearly and be influenced by a linear control action with probability \( 0\) in a generic observable space. This simple approximation of the soft robot is expected to fail in predicting the actual nonlinear dynamics of the soft robot. In fact, the authors in [27] observe this in their results and try to remedy their model by accounting for an additional term.

That additional term is called the Koopman static operator. Using the static inputs computed to learn the static SSM dictionaries, we learn this static operator by fitting a linear mapping \( \mathbf{G}\) between steady states in workspace \( \mathbf{z}_{s}\) and static inputs \( \bar{\mathbf{u}}\) , this is defined as:

They reformulate the optimization problem eq.(1) as a linear quadratic regulator (LQR) with a linear time-invariant dynamical system constraint. They then treat the finite-horizon problem as an infinite-horizon problem and find exact solutions for the desired control inputs required to track the target \( \boldsymbol{\Gamma}(t)\) by solving an algebraic Riccati-difference equation for the pregain matrix \( \mathbf{K}\) .

Since this resulted in inaccurate results, they added a static correction to the above as,

Once this term was added, [27] reported accurate tracking, and their results also indicated that the dynamic part had negligible contributions. Their method can only be reasoned to work when the target trajectory they wish to control is slowly varying and the robot’s dynamics is purely linear. They are invariably using these assumptions to first move the robot along slowly varying circular trajectories and then speed up the robot once they have already latched onto the target during the slow run. These methods will ultimately fail to achieve inertial dynamics if from the very start one challenged the Koopman pregain LQR controller to track a fastly accelerated target.

[27] also employ significant simplifications to the optimal control problem. Due to this, we expect their method to fail in tracking various targets precisely. This is exactly what we observed for the examples considered in this paper. Also, we note that the construct in [27] does not take into account the state and input constraints.

Generating the common training data for the aSSM and Koopman approaches took about an hour and a half for both of our examples. This time includes the numerical simulation of the systems involved, which would not be needed in an experimental implementation. Therefore, an hour and a half is an upper estimate of what the data collection time would be in an experimental setting. Once the training data were available, the model was trained in about 25 seconds for the aSSM approach and about 1 second for the Koopman-based approach. These training times will be the same in an experimental implementation. For TPWL, generating the training data took less than an hour, whereas training the model took about thirty seconds. None of these times is extensive, but the results from the different methods show large variations. Specifically, the accuracy of the aSSM-based models is consistently orders of magnitude higher than that of the Koopman and TPWL models, even though the latter two are higher-dimensional. All this is not surprising, given that aSSM-reduction targets a specific invariant manifold family that is known to exist in these highly nonlinear systems. In contrast, the Koopman-based and TPWL approaches fit linear models to those systems.

We isolate \( 1000\) random static steady states from the collected static input training data. In Fig. 46a, we plot the static steady states in the short arm’s workspace. We observe that these states accumulate along a tapered sphere with maximum height \( 1 \text{ [m]}\) and radial span \( 0.8 \text{ [m]}\) . The tapering is due to under-actuation and limited torque magnitudes. For these reasons, the maximum lobe sections of our planned Trifolium target track cannot be readily reached by the short elastic arm. To assess closed-loop performance, we quantify acceptable error regions by plotting an operational tolerance tube around the target track of radius \( 15 \text{ [cm]}\) . Similar conclusions about the steady-state distributions hold for the long elastic arm, where the tapered sphere has a wider span than in the short arm.

Pan/Zoom

Pan/Zoom

Pan/Zoom

Pan/Zoom

We generate responses of the short and long soft arms to a randomly varying 30-second control input, with a slowness measure \( r_s = 0.75\) , using the undeformed arm configuration as the initial condition. We pass this control input and initial condition to the learned \( 6\) D aSSM-reduced model in the form (13) and output a prediction in the robot’s workspace. In Fig. 46b and 46c, we plot the prediction of the aSSM-reduced model and the true response of the short and long soft arms. Our \( 6\) D aSSM-reduced model accurately predicts the transient dynamics from the origin to the slow steady state, as well as the slow steady state itself. The prediction closely matches the ground truth with an NMTE \( \approx 9 \%\) for the short soft arm and an NMTE \( \approx 11 \%\) for the long soft arm.

We add noise sampled from a bounded Gaussian distribution with an upper bound of \( 10 \text{ [mm]}\) to the feedback term in the aSSM-reduced MPC scheme eq.(26). In Fig. 46d, we plot the closed-loop result for the 6D aSSM-reduced model in the presence of feedback noise. The closed-loop prediction has visible jitter but overall tracks the target and remains within the operational tolerance tube. Here, the operational tolerance tube is enlarged to a radius of 18.5 [cm] to account for the noise uncertainty. The normalized ISE increases by a small factor of 0.01 compared to the case without noise presented in Section 6.3.