A simultaneous space‐time discretization approach to the inverse dynamics of geometrically exact strings

A continuous space‐time Galerkin method is newly proposed for the numerical solution of inverse dynamics problems. The proposed space‐time finite element method is combined with servo‐constraints to partially prescribe the motion of the underlying mechanical system. The new approach to the feedforward control of infinite‐dimensional mechanical systems is motivated by the classical method of characteristics. In particular, it is shown that the simultaneous space‐time discretization is much better suited to solve the inverse dynamics problem than the semi‐discretization approach commonly applied in structural dynamics. Representative numerical examples dealing with elastic strings undergoing large deformations demonstrate the capabilities of the newly devised space‐time finite element method.


INTRODUCTION
The present work deals with the numerical solution of inverse dynamics problems associated with the feedforward control of geometrically exact strings. In this specific type of inverse dynamics problem the motion of the system is assumed to be partially prescribed and the goal is to determine the actuating forces along with the overall motion of the system such that the partially prescribed motion is realized. In the context of infinite-dimensional models for strings flatness-based methods have been proposed in References 1-3 to solve the inverse dynamics problem (see References 4,5 for more background on differentially flat systems). The flatness-based methods typically rely on mechanical modeling assumptions such as small deformations 1,3 or inextensibility. 2 Concerning the feedforward control of underactuated finite-dimensional (or discrete) systems, servo-constraints provide an alternative to flatness-based methods. 6,7 For example, servo constraints have been successfully applied to underactuated mechanical systems such as cranes [8][9][10] and manipulators with passive joints. [11][12][13] In servo-constraint problems the motion of the system is partially specified by appending servo-constraints to the underlying equations of motion. Accordingly, the motion of servo-constrained discrete systems is governed by differential-algebraic equations (DAEs). In comparison to systems subject to standard holonomic constraints, servo-constraints * often lead to an increased index of the corresponding DAEs, see References 14-17 for more background on the (differentiation) index of DAEs. In this case appropriate index reduction techniques need be applied prior to the numerical time integration. 8,18 F I G U R E 1 The planar model of an overhead trolley crane F I G U R E 2 Elastic string (left) as sub-system of a multibody system (right): Cooperative transportation of a rigid body through four elastic strings actuated by external forces f k 0 , k ∈ {1, … , 4}. In the inverse dynamics problem the motion of the rigid body is prescribed and the goal is to determine the motion of the whole multibody system along with the actuating forces f k 0 . The prescribed motion of the rigid body gives rise to the boundary conditions (trajectory and contact force f k L ) for each string at its point of attachment to the rigid body A prototypical example for the type of inverse dynamics problems under consideration is the overhead trolley crane depicted in Figure 1 (see References 8,18 for further details). The mechanical model for the crane has three degrees of freedom associated with the generalized coordinates (s, l, ). Since only two actuating generalized forces (F, M) are present, the discrete system at hand is underactuated. The goal of the inverse dynamics problem is to determine the control inputs (F, M) such that the load (point mass m) traces a prescribed trajectory . For that purpose typically rest-to-rest maneuvers are considered. To realize the prescribed motion of the load, servo-constraints can be appended to the equations of motion pertaining to the crane. The resulting servo-constraint problem is governed by DAEs with differentiation index equal to five. Since the inverse dynamics problem is known to be differentially flat, 19 the inverse problem is well-posed.
In the present work we extend the notion of servo-constraint problems to the realm of infinite-dimensional systems. In particular, we focus on the inverse dynamics of elastic strings undergoing finite deformations. This specific problem can be viewed as sub-problem of the cooperative transportation of a cable-suspended load (rigid body) by means of unmanned aerial vehicles 3,20 (Figure 2). In the inverse dynamics problem of the multibody system the motion of the rigid body is prescribed and gives rise to prescribed boundary values (trajectory and contact force f k L ) for each string associated with its attachment point on the rigid body. In essence, the feedforward control of the whole multibody system boils down to the associated inverse dynamics problem for each elastic string. 21 The specific type of inverse dynamics problem for the elastic string entails non-standard space-time boundary conditions: while both the trajectory and the contact force are prescribed at one boundary, the external force acting on the other F I G U R E 3 Illustration of the geometrically exact string boundary is to be determined. As mentioned above, this type of inverse dynamics problem has been previously addressed in the context of flatness-based methods. It is important to realize that differential flatness of the problem implies that the inverse dynamics problem at hand is well-posed.
In the present work we newly propose to combine the notion of servo-constraints with a continuous space-time Galerkin approach to solve the specific space-time boundary value problem (BVP) at hand. In contrast to the previously developed flatness-based methods, our approach does not necessitate any simplifying modeling assumptions and thus makes possible to employ (truly nonlinear) geometrically exact strings. 22 In contrast to the semi-discretization approach commonly applied in structural dynamics (see, e.g., References 23,24), the simultaneous space-time discretization approach proves to be the natural choice for the solution of the inverse dynamics problem at hand. A readable account of previous developments in the field of space-time finite elements can be found in the recent work, 25 see also  We further show that the classical method of characteristics (see, e.g., References 30,31) can also be applied to solve the present space-time BVP. In particular, the method of characteristics provides further insight into the present space-time BVP, explaining why the simultaneous space-time discretization is much better suited to solve the inverse problem than the standard semi-discretization approach.
An outline of the rest of the article is as follows. Section 2 contains a summary of the formulation of geometrically exact strings along with the description of the inverse dynamics problem to be solved. In Section 3 the common semi-discretization approach is addressed and shown to be unsuited for the inverse dynamics problem at hand. The simultaneous space-time discretization is dealt with in Section 4. In particular, the classical method of characteristics is dealt with first, and paves the way to the newly proposed space-time finite element approach. After the presentation of representative numerical investigations in Section 5, conclusions are drawn in Section 6.

Geometrically exact strings
We focus on elastic strings undergoing finite deformations. An exhaustive account of the underlying geometrically exact description of strings can be found in the book by Antman. 22 To identify material points lying on the reference curve of the string we make use of the arc-length in the reference configuration s ∈ S = [0, L] ⊂ R ( Figure 3). In this connection, we assume a stress-free reference configuration in which the string has length L.
The placement of material point s ∈ S at time t ∈ T = [0, ∞) is characterized by the deformation map r(s, t) ∶ S × T  → R d , where d ∈ {1, 2, 3} is the spatial dimension. In the sequel it proves convenient to introduce the space-time domain Balance of linear momentum gives rise to the underlying equations of motion which take the form of the following system of partial-differential equations (PDEs) for all (s, t) ∈ Ω. Here, n(s, t) ∶ Ω  → R d represents the contact force at (s, t) ∈ Ω, b(s, t) ∶ Ω  → R d denotes the body force per unit reference length at (s, t) ∈ Ω and ( A)(s) ∶ S  → R + is the mass density per unit reference length at s. In this connection, stands for the density per unit reference volume and A refers to the cross-sectional area of the string in the reference configuration. The body force per unit reference length is given by where g is the vector of gravitational acceleration.
Since the string has no bending stiffness, the contact forces in the string need be oriented tangentially to the string. That is, The tension N(s, t) at (s, t) ∈ Ω is determined by a frame-indifferent constitutive relation of the form where the stretch at (s, t) ∈ Ω is defined by The total actual length of the string follows from Note that the stretch represents the local ratio of the deformed to the reference length of the string. A physically meaningful elastic constitutive law needs to satisfy the restrictionsN( , s) → ∞ as → ∞ andN( , s) → −∞ as → 0. In addition to that, we assume a natural reference configuration, implyingN(1, s) = 0 (i.e., vanishing tension in the reference configuration).
Example 1 (Hyperelastic material model). In the present work we apply a hyperelastic constitutive model for a homogeneous rope based on the stored energy Û( ) = EA 4 ( 2 − 2 ln − 1). This stored energy represents a simple model for uniaxial rubber-type material behavior. The tension in the string follows fromN( ) = Û ′ ( ) leading to the constitutive Linearization of the last equation about the natural reference configuration yields the linearized constitutive relationship Thus EA can be regarded as axial stiffness in the reference configuration.

Inverse dynamics of nonlinear strings
We assume that the motion of the elastic string is partially specified at its right boundary at s = L ( Figure 4). That is, the placement of the free end of the string at s = L is specified for all time, leading to the boundary conditions F I G U R E 4 Illustration of the inverse dynamics problem at hand where (t) ∈ R d is a prescribed function of time. The goal is now to find the actuating force vector f 0 (t) ∈ R d acting on the left boundary at s = 0 such that the partially prescribed motion of the string is realized. The corresponding boundary condition can be written as The resulting feedforward control problem gives rise to an initial boundary value problem which can be stated as follows: Here, the two abbreviations have been introduced. The PDEs (11) 1 result from substituting (3) into the equations of motion (1). Note that these PDEs can be classified as second-order quasilinear hyperbolic system in one dimension for r(s, t) (see Antman 22 ). The nonlinearity of (11) 1 is reflected in the dependence of B(s, t) on the unknown vector r(s, t), compare (12) 2 together with (4) and (5). Moreover, (11) 2 correspond to the boundary conditions (10) and (9), respectively. Accordingly, if the right end of the string is free, we have f L (t) = 0. Note, however, that the force vector f L (t) has been introduced to take into account more general scenarios such as those mentioned in Remarks 1 and 2.
Initial conditions are account for by (11) 3 , where f 0 ∈ R d and r 0 , v 0 ∈ R d are prescribed. Feedforward control problems often deal with rest-to-rest maneuvers. There, the solution of the corresponding equilibrium problem supplies the initial values f 0 , r 0 and v 0 = 0. In the more general case in which the initial configuration of the string is not at rest, the initial values for the inverse problem can be obtained by solving the corresponding forward dynamics problem.
Remark 1 (Additional point mass). An additional point mass M might be attached to the right boundary of the string (at s = L). In this case, the last boundary condition in (11) 2 would have to be replaced by where 2 t r(L, t) =̈(t) and, as before, g denotes the vector of gravitational acceleration.

F I G U R E 5 The linear-elastic bar
Remark 2 (Multibody system). If the elastic string is a sub-system belonging to the inverse dynamics problem of a multibody system such as that shown in Figure 2, the force vector f L (t) in the last boundary condition in (11) 2 is a prescribed function of time (see Reference 21 for further details).
Example 2 (Linear-elastic bar). Concerning small longitudinal deformations about a straight reference configuration, we recover the model of a linear-elastic bar ( Figure 5). Accordingly, consider the one-parameter family of configurations r (s, t) = r (s, t)e, where e is a unit vector and r is a power series expansion of the form Here, r 0 = s specifies the placement of material points in the straight reference configuration and d d | | | =0 r (s, t) = u(s, t) characterizes small longitudinal displacements. Moreover, the axial strains can be obtained via Similarly, the linearized stretch yields = 1 + s u, so that the linearized constitutive relation (8) can be written asN lin = EA . Eventually, linearizing the equations of motion (11) 1 of the string about the straight reference configuration and neglecting the body forces yields Here, for the sake of clarity and without loss of generality, a free end at s = L is assumed. In (14) 1 , c stands for the constant wave propagation speed defined by There exists an analytical solution to problem (14) that will serve as reference for the numerical methods developed in the sequel.
Remark 3 (Analytical solution). An analytical solution to problem (14) is based on d'Alembert's formula (see, e.g., Reference 32 or 33 (ch. 16)). It can be easily verified that satisfies PDE (14) 1 . Now, the first boundary condition in (14) 2 yields f (t) = EA s u(0, t), where (16) implies Inserting (16) into the second boundary condition, u(L, t) = (t), and subsequently differentiating with respect to time yields Similarly, the third boundary condition, s u(L, t) = 0, gives Adding (18) and (19) Evaluating the last two equations at t = t − c −1 L and t = t + c −1 L, respectively, leads to Substituting from (20) into (17), the first boundary condition eventually yields the result providing the actuating force at the left boundary of the bar in terms of the prescribed displacement of the right boundary of the bar.

SEQUENTIAL SPACE-TIME DISCRETIZATION
In nonlinear structural dynamics the solution of initial BVPs is commonly performed by applying a sequential space-time discretization (see, e.g., References 23,24). The so-called method of lines is typically based on a discretization in space by means of finite elements, followed by a subsequent discretization in time by means of finite differences. Following this common procedure we show that the inverse dynamics problem under consideration can be transferred to discrete equations of motion subjected to servo constraints. For related discrete servo-constraint problems we refer to References 7,8,34.

Semi-discrete equations of motion
In this section we provide a short outline of the common space discretization of initial BVPs by applying the finite element method. For that purpose we consider the pure Neumann problem † associated with the initial BVPs (11). Multiplying (11) 1 by a sufficiently smooth test function w(s) ∈ R d and integrating subsequently over the spatial domain S yields Applying integration by parts to the second integral on the left-hand side and taking into account the Neumann boundary conditions leads to the weak form The weak form can now be discretized in space by applying finite element approximations to the vector-valued test function w(s) and the trial function r(s, t). We confine our attention to piecewise linear approximations based on Lagrangian shape functions (cf. Reference 35). The corresponding interpolations read

F I G U R E 6 Discretized spatial domain
where L i (s) are linear Lagrangian shape functions related to the mesh sketched in Figure 6. The associated nodal values Inserting the finite element approximations into weak form (23) and assuming for simplicity f L (t) = 0 yields the semi-discrete equations of motion Here, q(t) is the nodal configuration vector that contains the nodal position vectors at time t ∈ T of the discrete problem at hand. That is, Furthermore, the nodal contributions to mass matrix M, internal force vector f int (q), and external force vector f ext (t) are given by, respectively, where I d is the d × d identity matrix and b(s, t) is given by (2). Matrix B 0 in (25) is of Boolean type and essentially links force f 0 (t) to the first node lying at the left boundary of the string, accounting for the virtual work contribution q ⋅ B T 0 f 0 (t) = w 1 ⋅ f 0 (t) emerging from the right-hand side of weak form (23). Here, q contains the nodal values w i in analogy to (26).
The semi-discrete equations of motion (25) constitute a system of nonlinear ordinary differential equations (ODEs) of second-order. In the standard forward dynamics problem the forcing terms f ext (t) and f 0 (t) are prescribed functions of time and the ODEs (25) can be solved by applying common time-stepping schemes.
Remark 4 (Additional point mass). If an additional point mass M is attached to the right end of the string (Remark 1), has to be taken into account in weak form (23). Correspondingly, in the semi-discrete formulation, the following entries of mass matrix (27) 1

Servo-constraints
In the inverse dynamics problem stated in Section 2.2, the actuating force f 0 (t) plays the role of an unknown which needs be determined by accounting for the partially specified motion of the string at its right boundary. This can be accomplished by appending servo constraints of the form r p+1 (t) = (t) to the semi-discrete equations of motion (25). The servo-constraints can be recast in the form F I G U R E 7 Sketch of the semi-discrete string model subject to servo-constraints where C is a Boolean matrix extracting from the nodal configuration vector q(t) the entries associated with node p + 1 lying at the right boundary of the discrete string ( Figure 6). A sketch of the semi-discrete string model subject to servo-constraints is depicted in Figure 7.
The discrete servo-constraint problem is governed by the differential equations (25) together with the algebraic equations (28). Accordingly, the servo-constraint problem gives rise to differential-algebraic equations (DAEs). However, these DAEs have to be distinguished from the DAEs governing the motion of discrete mechanical systems subject to holonomic constraints. If the constraints were holonomic, the term B T 0 f 0 in (25) would correspond to constraint forces with associated Lagrange multiplier f 0 for the enforcement of holonomic constraints ‡ B 0 q(t) = (t). Moreover, in the standard case of holonomic systems the DAEs have (differentiation) index d = 3 and are thus amenable to a direct numerical solution (see, e.g., Reference 37).
Often the DAEs related to servo-constraint problems have index d ≤ 5 7 and can still be integrated numerically by applying appropriate index reduction techniques. 8,18 In contrast, however, the index of the present DAEs turns out to increase with the number of finite elements used for the space discretization. This is illustrated next with the example of the linear-elastic bar.

Example 3.
In the case of the linear-elastic bar (Example 2), the semi-discrete equations of motion (25) in conjunction with the servo-constraints (28) boil down to the DAEs For simplicity, we consider a lumped mass matrix with entries M ij = M i ij . Considering the last two nodes p and p + 1 associated with the p th element ( Figure 6), (29) yield Here, the summation convention does not apply. Since (t) is prescribed, it can be deduced from the last equation that the displacement of node p is given by Consequently, where (n) = d n ∕dt n and function a p follows directly from differentiating (31) twice with respect to time. Repeating this procedure for the elements p − 1, … , p − n, reveals that Setting n = p − 1 leads to so that for i = 1, (29) 1 yields One more time differentiation of the last equation reveals that the time derivative of the algebraic variable f pertaining to the DAEs (29) can be obtained after 2p + 3 time derivatives of the algebraic constraint (29) 2 . This implies that the DAEs (29) have differentiation index d = 2p + 3 (cf. References 16,17). It can be shown that this result also holds true in case of a consistent mass matrix. Accordingly, the index increases with the number p of finite elements used for the space discretization of the bar.
Remark 5 (Differential flatness). Since all of the unknowns of the DAEs (29) can be expressed in terms of and derivatives thereof, the present problem enjoys the property of differential flatness. 4,5 In particular, the nodal displacement u p+1 plays the role of a flat output. Analogous results have been obtained for the articulated mass point systems dealt with in References 6,34,38. Note that the simple result for the flat output in Example 3hinges on the mass matrix being diagonal. The situation gets more intricate in case of a consistent mass matrix or higher-order finite elements.
To summarize, the sequential space-time discretization approach to solve the inverse dynamics problem at hand yields DAEs whose index tends to be excessively high thus hindering a stable numerical solution. In addition to that, the demands on the smoothness of the prescribed trajectory tend to be excessively high as well.

SIMULTANEOUS SPACE-TIME DISCRETIZATION
To circumvent the problems resulting from the sequential space-time discretization outlined in Section 3, we propose a simultaneous space-time discretization for the solution of the inverse dynamics problem at hand. Our approach is guided by the classical method of characteristics which is known to be applicable to solve the wave equation.

Method of characteristics
In this section we apply the method of characteristics to solve the inverse dynamics problem (11). The method of characteristics is based on a geometric interpretation of quasi-linear partial differential equations (see References 30,31,39-41 for more background information).
In a first step we convert the second-order PDEs (11) 1 into a system of first-order PDEs. For that purpose we introduce the velocity space-time field v(s, t) = t r(s, t) along with the tangent vector p(s, t) = s r(s, t). Now, (11) 1 can be rewritten as where in the last equation the equality of mixed derivatives, 2 ts r = 2 st r, has been used. We aim at employing the state vector To this end we make use of the relationship Taking into account that B(s, t) =B(r(s, t), s), the last term on the right-hand side of (38) can be rewritten as Substituting from (38) and (39) into (36) where Equation (36) 1 together with (40) give rise to the following system of first-order PDEs where D = Assume there exists a curve s = k(t) along which a solutionẑ(t) = z(k(t), t) is given. Then this curve is called a characteristic if the partial derivatives of the solution cannot be uniquely determined through information along this curve. Accordingly, inserting Since, by assumption, s z can not be determined by knowledge of z on the curve, the condition has to be satisfied. Consequently, k ′ (t) coincide with the eigenvalues of E relative to D (i.e., the eigenvalues of D −1 E). Correspondingly, there are 2d characteristics associated with the eigenvalues i (i = 1, … , 2d). In particular, for the hyperelastic material model introduced in Example 1, matrix (41) assumes the form so that condition (45) yields the pairwise solution for , and k ′ 2 (t) = ±c for d = 2 and for d = 3. In the above equations c denotes the constant wave propagation speed introduced in (15).
In the sequel we illustrate the application of the method of characteristics to the inverse dynamics problem with Examples 4 (linear case) and 5 (nonlinear case). Our approach can be traced back to the graphical-numerical solution method due to Massau (see Reference 42).
Remark 6. Apart from the one-dimensional case (d = 1), the eigenvalues (48) and (49) are not all real-valued anymore for < 1. This conforms with the fact that strings can sustain tensile forces but not compressive ones.
Example 4 (Linear-elastic bar). For the linear-elastic bar in one dimension (d = 1) introduced in Example 2, the second-order PDE (14) 1 is replaced by the system of first-order PDEs (42) based on the state vector The characteristics are characterized by the linear version of (47) given by Accordingly, the slope of the characteristic lines in (s, t) space is determined by the constant wave propagation speed introduced in (15). In addition to condition (45), problem (44) gives rise to two further conditions Here,v(t) = v(k(t), t) andn(t) = n(k(t), t) stand for, respectively, the velocity and the normal force along the characteristic lines s = k(t) whose slope k ′ (t) is given by (51). The two conditions in (52) eventually yield one independent equation of the form The ODEs (53) can be solved along the characteristic lines specified by (51). Altogether, application of the method of characteristics converts the scalar second-order PDE (14) 1 into a system of four ODEs given by (51) and (53). The ODEs are supplemented with boundary conditions and initial conditions v(s, 0) = v 0 and n(s, 0) = 0. For the linear problem at hand the ODEs (51) and (53) are decoupled. Accordingly, (51) can be used in a first step to set up the characteristic net consisting of straight characteristic lines specified by (51). A sample characteristic net thus obtained is shown in Figure 8. To each node (s Q , t Q ) ∈ Ω of the characteristic net there is associated a nodal state vector z Q = (v Q , n Q ). The unknown nodal states can be determined by integrating (53) along the characteristic lines. This procedure is illustrated with Figure 9 in which three nodes Q and P j , j ∈ {1, 2} are considered. The three nodes correspond to vertices on the characteristic net and thus comply with (51) in the sense that Characteristic net along with a representative node Q located at (s Q , t Q ) ∈ Ω F I G U R E 9 Three nodes Q, P 1 , and P 2 , located on the characteristic net with associated states (v Q , n Q ), (v P 1 , n P 1 ), and (v P 2 , n P 2 ). The characteristic lines C 1 and C 2 are associated with the slopes ds∕dt = c and ds∕dt = −c, respectively is satisfied for j ∈ {1, 2}. The numerical integration of (53) relies on the finite difference scheme for j ∈ {1, 2}. Repeatedly applying this procedure to the characteristic net in Figure 8 Taking the time derivative of (55) along the characteristic lines yields where (53) has been used, which is valid along the respective characteristic line. Accordingly, the functions (55) ( = 1, 2) are indeed constant along the characteristics associated with eigenvalues .
Example 5 (Nonlinear string in 2d). For the nonlinear string in two dimensions (d = 2), the second-order PDE (11) 1 is replaced by the system of first-order PDEs (42) based on the matrices in (43) and (46). The state vector (37) is comprised of the 2-vectors v and n. Two pairs of characteristics can be distinguished which are characterized by the slopes in (48). Correspondingly, two (i = 1, 2) pairs of characteristics are defined by where, with regard to (48), Each pair of characteristics has j ∈ {1, 2} members, where j = 1 refers to the forward propagating direction and j = 2 to the backward propagating direction. In addition to (57), (44) gives rise to the following ODEs along the characteristics wherev = v(k(t), t) andn = n(k(t), t) denote, respectively, the velocity and the contact force along the respective characteristic curve. Moreover, Equations (57) and (59) constitute a system of coupled nonlinear ODEs which can be discretized as follows: for i, j ∈ {1, 2}. In analogy to Example 4, vertices of the characteristic net are denoted by (s A , t A ), while the corresponding nodal states are denoted by (v A , n A ). In contrast to Example 4, the characteristic net can not be determined beforehand, because c i in (60) 1 depends on the state, see (58). Consequently, the location of the nodal net points, (s Q , t Q ) and (s P ij , t P ij ), in general needs be determined through the solution procedure, as well as the nodal states. Scheme (60) is further illustrated with Figure 10, in which the points Q and P ij (i, j ∈ {1, 2}) located on the respective characteristics are displayed. This procedure can be repeatedly applied within the space-time domain of interest ( Figure 11). After assembly, a coupled system of nonlinear algebraic equations is obtained which can be solved iteratively F I G U R E 10 Integration along the characteristics C ij (i, j ∈ {1, 2}) according to (60). The characteristics C ij are associated with the slopes (ds∕dt) ij in (57)

F I G U R E 11
Characteristic net in the space-time domain by means of Newton's method. Prior to the solution the boundary and initial conditions emanating from (11) 2,3 need be taken into account.
Remark 9 (Pre-and post-actuation phases). Note that similar to the linear setting (Remark 7), the solution of the inverse dynamics problem requires pre-and post-actuation phases. This is again indicated with the shaded area in Figure 11. The parameters t 0 , t f , and T have the same meaning as in Remark 7.

Space-time finite element approach
The method of characteristics dealt with in Section 4.1 can be viewed as simultaneous space-time discretization approach to the inverse dynamics problem (11). Aiming at an alternative, more systematic, approach we next propose a simultaneous space-time finite element method to solve the inverse dynamics problem at hand. Introducing the velocity field v(s, t) in space-time, the PDEs (11) 1 can be recast in the form Scalar multiplying the equations in (61) by sufficiently smooth test functions w 1 , w 2 ∶ Ω  → R d and subsequently integrating over the space-time domain Ω = S × T yields the following integral formulation: Integrating by parts the second term on the left-hand side of (62) 2 and taking into account the boundary conditions (11) 2 yields where boundaries Γ 0 = {0} × T and Γ L = {L} × T of space-time domain Ω = S × T have been introduced (see Figure 12). Concerning the prescribed motion of the right end of the string (at s = L), we make use of weakly enforced servo-constraints of the form where w 3 ∶ Γ L  → R d is a third test function. Now the weak formulation of the inverse dynamics problem at hand can be stated as follows: Given b, f 0 , r 0 , and v 0 , find r ∈ V 1 , v ∈ V 2 , and f 0 ∈ V 3 , such that for arbitrary test functions w ∈ W ( = 1, 2) and w 3 ∈ W 3 . The prescribed initial conditions (11) 3 give rise to the Dirichlet-type boundary Ω 0 of the space-time domain Ω ( Figure 12). Accordingly, we have where n nodes is the total number of nodes given rise to the index set  nodes = {1, … , n nodes }. Moreover, index set  0 ⊂  nodes contains the nodes lying on Γ 0 . Similarly,  L ⊂  nodes contains the nodes lying on Γ L . The test functions are approximated by

Resulting algebraic formulation
In the present work we confine our attention to standard low-order isoparametric finite elements. In particular, L i denote linear Lagrangian shape functions also employed in Section 3.1, while N a denote bilinear Lagrangian shape functions. Inserting the above finite element approximations into weak form (65) yields the algebraic system of equations where In (66), the unknown nodal quantities r a , v a , and f 0 i have been assembled in corresponding nodal system vectors Q, V, and  0 . Elimination of the nodal velocity vector V yields the residual vector (68) The resulting algebraic system of nonlinear equations, (Q,  0 ) = 0, can be solved iteratively for the nodal unknowns (Q,  0 ) by applying Newton's method. In this connection, the iteration matrix is given by where DF int (Q) results from (67) 3 together with the constitutive term (12) 2 . A straightforward calculation taking into account the material model outlined in Example 1 yields the nodal contribution to DF int (Q) given by where matrix H has been defined in (41). Note that p in (41) has now to be substituted with r h .
Remark 10 (Additional point mass). If an additional point mass M is attached to the right end of the string (Remark 1) , has to be taken into account in weak form (65) 2 . Correspondingly, in the discrete formulation, the following entries of matrices (67) 1 and (67) 5 need be modified according to

Recursive implementation
The newly devised space-time finite element method yields the system of algebraic equations (66) which can be solved by applying the iterative procedure outlined in Section 4.2.1. This implies the simultaneous solution of all the unknowns resulting from the space-time discretization. We next show that the solution can alternatively be attained by applying a recursive scheme which relies on the decomposition of Ω into N time-space slabs Ω n = (s n−1 , s n ) × T, n = 1, … , N ( Figure 13). Focusing on one representative slab Ω n , the recursive solution procedure relies on the weak form where Γ n = {s n } × T. Note that Γ N = Γ L . In essence, weak form (71) emanates from the original formulation (65) by restricting the space-time domain to Ω n . In this connection, (71) 3,4 have been introduced to facilitate a convenient description of the recursive solution procedure. The overall space-time approximation remains the same as before. The approximation of the newly introduced intermediate quantities f n (n = 1, … , N − 1) is given by where index set  n contains the nodes lying on Γ n . In this way, and f n (n = 1, … , N − 1) have to be prescribed in the same way as f 0 , see Section 2.2. Similarly, we introduce It can be easily shown that the recursive application of (71) for n = 1, … , N is equivalent to the original space-time method, provided that the approximation of n in (71) 3,4 conforms with the continuity of the displacement approximation.
To ensure this, we choose so that (71) 3,4 enforce r h| |Γ n = h n and r h| |Γ n−1 = h n−1 . Now, the equivalence between (71) 1 and (65) 1 follows from the property Similarly, the equivalence between (71) 2 and (65) 2 follows by additionally taking into account Moreover, (71) 3 contains (65) 3 for n = N under the provision that the prescribed trajectory on Γ L is based on nodal interpolation of the type (73). Now the proposed space-stepping algorithm results from the recursive application of the discretized weak form (71) in backward space direction, that is, for n = N, N − 1, … , 1 as shown in Table 1.
Note that for n = N the given data follows from the prescribed data on Γ L leading to h N = h and f h N = f h L , where the prescribed data and f L is assumed to be given in interpolated form of the type (73) and (72), respectively.
The implementation of the recursive scheme is again based on the algebraic formulation outlined in Section 4.2.1. However, in each step of the recursive scheme all of the algebraic quantities are now confined to the respective slab Ω n (see also Remark 12 below).
Accordingly, f n = n| Γ n and f n−1 = n| Γ n−1 , so that f n and f n−1 correspond to the contact forces in the string along Γ n and Γ n−1 , respectively.
Remark 12. (Numerical effort). Assuming a regular discretization of the space-time domain Ω = S × T relying on n s nodes in space direction and n t nodes in time direction amounts to a total number of nodes n nodes = n s ⋅ n t . The simultaneous solution for all unknowns outlined in Section 4.2.1 relies on the iteration matrix (69), whose dimension is equal to (n s + 1) ⋅ (n t − 1) ⋅ d. In contrast to that, the dimension of the iteration matrix corresponding to each of the N = n s − 1 steps of the recursive solution procedure is equal to 3 ⋅ (n t − 1) ⋅ d.
Remark 13. The semi-discrete approach described in Section 3 can be linked to the space-time finite element method developed above. To this end, we consider a Galerkin-based discretization in time of the semi-discrete formulation. In particular, we introduce the following approximations of the nodal position and velocity vectors For conciseness the summation convention applies throughout this remark. Moreover, for simplicity we again employ Lagrangian shape functions in (74). In addition to that, we introduce nodal weighting functions of the form where M A are basis functions to be specified below. Now the semi-discrete equations (25) together with the servo-constraints (28) can be recast in the weighted Galerkin form Note that (76) 1 serves the purpose of introducing the nodal velocities v h i . Time-stepping schemes typically applied in the context of the semi-discrete formulation can be recovered from (76) by choosing specific discontinuous shape functions M A in (75), along with associated quadrature formulas for the evaluation of the time integrals. This approach gives rise to the so-called continuous Galerkin method § (see Reference 33). For example, focusing on the first two equations in (76) emanating from the underlying system of first-order ODEs, choosing constant shape functions M A and linear Lagrangian shape functions L j , along with the mid-point quadrature for the evaluation of the time integrals yields the mid-point rule (see, e.g., Reference 44).
In sharp contrast to that, the newly developed space-time finite element method can be recovered from (76) by choosing continuous shape functions M A . In particular, using again Lagrangian shape functions, that is, M A (t) = L A (t), (76) leads to the algebraic formulation (66), provided that • the shape functions N a (s, t) used in the above space-time finite element method have tensor product form (cf. References 45,46). That is, Referring to the local node numbering of a bi-linear master element (Figure 14), the node numbers in (77) are related to each other as depicted in Table 2.
In this connection, the one-dimensional linear Lagrangian shape functions are given by where 1 and 2 refer to the nodal coordinates, while stands for either s ∈ S or t ∈ T.
• rectangular shaped space-time finite elements are used leading to a structured mesh in space-time. Note that this precludes unstructured space-time meshes which are feasible with the space-time finite element method proposed above.
To further illustrate the connection between (76) and (66), we consider the first term in (76) 2 , which yields Here, M ij introduced in (27) has been used. Moreover, the relation between the nodal shape functions in Table 2 has been taken into account. The remaining terms in (76) can be linked to the corresponding terms in (66) in a similar manner.

NUMERICAL INVESTIGATIONS
In this section the inverse dynamics problem at hand is numerically solved by applying the newly developed space-time discretization methods. In particular, the method of characteristics, termed "MOC" (Section 4.1), and the space-time finite element method, termed "ST-FEM" (Section 4.2), are applied. ST-FEM relies on bi-linear isoparametric space-time elements. The evaluation of the integrals in (67) relies on standard Gaussian quadrature rules.

Linear-elastic bar
We start with the feedforward control of the linear-elastic bar. This problem makes possible the comparison of the numerical results delivered by both MOC and ST-FEM with the analytical solution (Remark 3). The prescribed trajectory of the right end of the bar is assumed to be given by for t 0 = 1 and t f = 3. Moreover, the remaining data for this problem are EA = 1, A = 1, L = 1.
In Figures 15 and 16 the numerical results for the actuating force f 0 (t) are compared to the reference solution obtain by applying formula (21). It can be observed that the numerical results of the two alternative schemes under consideration closely match the reference solution. The results of MOC rely on a total of 356 unknowns. Concerning ST-FEM, two alternative meshes comprised of rectangular bilinear elements have been applied. The first one consists of 5 × 15 = 75 elements, leading to a total number of 195 unknowns. The second one consists of 10 × 50 = 500 elements, leading to a total number of 1150 unknowns. Figure 17 depicts the numerical solution for the functions v(s, t) = t r(s, t) and p(s, t) = s r(s, t), respectively, computed with the method of characteristics. The corresponding characteristic net is also partially visualized in the s, t-plane. Figure 18 contains plots of the two functions h (s, t) introduced in (55) associated with the Riemann invariants (Remark 8). In particular, function h 1 (s, t) corresponds to the characteristic line C 1 whose slope is given by ds∕dt = k ′ (t) = c. Similarly, function h 2 (s, t) corresponds to the characteristic line C 2 whose slope is given by ds∕dt = k ′ (t) = −c. It can be observed that both functions are indeed constant along the respective characteristic line.
The numerical solution for the displacement field u(s, t) and the velocity field v(s, t) computed with the space-time finite element method are depicted in Figure 19. There, the corresponding finite element mesh is also partially visible in the s, t-plane.

Nonlinear elastic string in 2d
The next example deals with the planar motion of a geometrically exact string. In particular, the trajectory of the right end of the string is prescribed as straight line in the x, y-plane via where function (t) is given by (79) with t 0 = 2 and t f = 4. The remaining data is given by EA = 1, A = 1, L = 1, g = 9.81.
The initial values for the inverse dynamics problem (11) have been obtained by solving the equilibrium problem of the string by means of the semi-discrete model described in Section 3.1. For that purpose, 50 finite elements have been employed. The thus obtained vertical equilibrium configuration of the string is characterized by a total length (6) of l(0) = 9.985 and a suspension force of f 0 = −g e y . Figure  To get an impression of the overall motion of the string, snapshots in the x, y-plane and in space-time are shown in Figures 22 and 23, respectively. In addition to that, the stretch distribution over the space-time domain is depicted in Figure 24.
In the second part of this example an additional point mass M = 1 is attached to the right end of the string (cf. Remarks 1,4, and 10). The prescribed motion of the point mass is again governed by (80), with t 0 = 1 and t f = 3. The remaining data is specified by EA = 10, A = 1, L = 1, g = 9.81.
In a first step the equilibrium configuration of the string-mass system has been computed by applying 50 finite elements within the semi-discrete formulation (Section 3.1). Accordingly, the vertical equilibrium configuration providing the initial data for the subsequent inverse dynamics simulation is characterized by a total length (6) of l(0) = 3.258 and a suspension force of f 0 = −2g e y . The resulting components of the actuating force required to realize the partially specified motion are depicted in Figure 25. It can be seen that the results of both methods under investigation match each other closely. The numerical results of the space-time finite element method are based on a mesh containing 50 × 200 = 10,000 elements leading to a total number of 20,800 unknowns. Merely 4 Newton iterations are required to attain the converged solution. The corresponding results of the method of characteristics are based on a total number of 10,515 unknowns. The required number of Newton iterations amounts to 7. Again the data pertaining to the initial equilibrium configuration of the system has been used to initialize Newton's method throughout space-time. The overall motion of the string is illustrated with snapshots of successive configurations of the system in Figure 26. The total length (6) of the string versus time is depicted in Figure 27. In addition to that, the stretch distribution over the space-time domain is depicted in Figure 28.

Nonlinear elastic string in 3d
The final example deals with the spatial motion of the string. In particular, the point mass attached at the right end of the string (cf. Remarks 1,4, and 10) is supposed to move on a helix. To this end, the helix-shaped trajectory is prescribed by In Figure 29 the prescribed coordinates are plotted versus time. The remaining data for this example is specified by EA = 10, A = 1, L = 1, g = 9.81, M = 1, and z f = 5. In a first step the vertical equilibrium configuration of the system ( Figure 30) subject to gravity load is computed by applying the semi-discrete model of the string (cf. Section 3.1) consisting of 15 finite elements. Accordingly, the vertical equilibrium configuration is characterized by a total length (6) of l(0) = 3.258 and a suspension force of f 0 = −2g e z . The thus obtained equilibrium configuration provides the initial data for the inverse dynamics problem (11) to be solved subsequently. In this connection, the data is also used to initialize Newton's method through the space-time domain. The computed three components of the actuating force required to realize the partially prescribed motion of the system are depicted in Figure 31. Figure 32 displays the total length of the string versus time. The resulting motion of the system is illustrated with Figures 33 and 34. These results have been obtained by employing 15 × 149 = 2235 rectangular space-time finite elements amounting to a total number of 7599 unknowns. Merely 4 Newton iterations were required to reach the numerical solution.

5.3.1
Recursive solution procedure We next apply the recursive solution procedure described in Section 4.2.2. The data is kept the same as before. The number of finite elements in space direction is reduced from 15 to 10 in order to simplify the graphical illustration of the results. Accordingly, the space-time finite element mesh is comprised of 10 × 149 rectangular elements leading to N = 10 time-space slabs used in the recursive solution procedure. As expected, the numerical results are practically indistinguishable from those obtained by applying the simultaneous solution procedure (Section 4.2.1). In Figure 35 the contact forces and trajectories along the boundaries Γ n of the time-space slabs are shown. While the dimension of the iteration matrix pertaining to the simultaneous solution is equal to 5364, the dimension of the iteration matrix of the recursive solution procedure is equal to 1341 (cf. Remark 12). On average 3.2 Newton iterations were required in each of the 10 steps of the recursive scheme. In comparison, 4 Newton iterations were required to reach the solution of the simultaneous solution scheme.
Remark 14 (Perturbed data of the prescribed trajectory). According to Hadamard a problem is well-posed if there exists a solution to the problem which is unique and stable (see Reference 47 or 48 (definition 1.13)). Since existence and uniqueness of the solution to the present space-time boundary value problem are assured by the differential flatness of the system, stability still has to be ensured. The solution is called stable if it depends continuously on the given data-small perturbations of the prescribed servo-constraints on Γ L must not lead to unbounded results for the actuating force f 0 (t) on Γ 0 . This can be verified numerically by adding random noise to the given data of the trajectory such that where denotes the perturbed trajectory. The perturbed components of are shown in Figure 36. For comparison, the original components of given by (81) are depicted in Figure 29. Figure 37 shows the results for the actuating forces resulting from the solution of the inverse problem with the perturbed data. Comparison with Figure 31 indicates that the perturbed problem still yields results close to the original ones.

F I G U R E 37
Components of the actuating force obtained as solution of the inverse problem with noisy data for the prescribed trajectory at Γ L with ∶ t  → [0, ] and ∈ {0.02, 0.05}

CONCLUSIONS
In the present work a new approach to the numerical simulation of inverse dynamics problems dealing with the feedforward control of geometrically exact strings has been developed. The newly proposed approach consists of the simultaneous discretization in space and time of the underlying space-time BVP. In particular, two alternative methods have been presented. While the first approach is based on the classical method of characteristics, the second approach relies on a Galerkin-based space-time discretization of the underlying PDEs. In essence, the space-time finite element method has been guided by the method of characteristics in that it motivates the simultaneous discretization in space and time. Moreover, in the finite element approach the trajectory tracking condition has been imposed weakly by means of servo constraints. Previously, servo-constraints were employed successfully in the simulation approach to the feedforward control of underactuated discrete mechanical systems such as cranes and manipulators with passive joints (see, e.g., References 9,10,34). It is tempting to follow similar lines by applying finite elements in the common way to convert the continuous model of the string into an underactuated discrete system subject to servo-constraints. However, we have seen in Section 3 that the resulting system of DAEs is prone to an excessively high index which hinders the stable numerical integration.
Thus the successive space-time discretization approach commonly applied in structural dynamics turns out to be rather unsuited to solve the inverse dynamics problem at hand.
A characteristic feature of the present space-time Galerkin finite element approach are continuous test functions. This is in sharp contrast to previously developed space-time finite element methods that rely on discontinuous test functions thus making possible to divide the space-time domain into a sequence of space-time slabs eventually recovering the time-marching format commonly applied in structural dynamics (see Reference 49, (sec. 6.3) and references cited therein).
Of course, the simultaneous space-time discretization approach yields a comparatively large system of algebraic equations whose unknowns refer to the whole space-time domain pertaining to the inverse dynamics problem to be solved. This feature is shared with the alternative approach based on the method of characteristics. On the other hand, the numerical investigations have shown that the total number of Newton iterations required to yield a converged solution is surprisingly low, although the system is undergoing quite large deformations. This puts the relatively large number of unknowns into perspective.
The numerical effort can still be reduced considerably by dividing the space-time domain into a sequence of time-space slabs making possible a recursive solution of the inverse problem as has been shown in Section 4.2.2. In essence, the recursive solution procedure does not change the structure of the algebraic system to be solved in each step but significantly reduces its size. Indeed, we have shown that both the recursive solution procedure and the simultaneous one yield equivalent results. While the recursive scheme leads to a reduction of the numerical effort, the simultaneous scheme is more general in that it makes possible to use unstructured space-time meshes.
If the actuating forces obtained from the inverse dynamics simulation are fed into a forward dynamics simulation based on the corresponding semi-discrete system in conjunction with a standard time-stepping scheme such as the mid-point rule, the trajectory previously prescribed in the inverse dynamics problem is accurately tracked. Applying this procedure to the 3d problem dealt with in Section 5.3, employing 149 time steps leads to an average of 2.85 Newton iterations per time step. Thus a total number of 425 Newton iterations is required in the forward dynamics simulation which is to be contrasted with merely 4 Newton iterations required to solve the inverse dynamics problem.
Not surprisingly, the method of characteristics is better adapted to the inverse dynamics problem than the space-time finite element method in the sense that it generally requires less unknowns to yield comparatively accurate results. In addition to that, the method of characteristics provides valuable insights into the inverse dynamics problem at hand. In particular, it explains why pre-and post-actuating phases are required to solve the inverse problem (cf. Remarks 7 and 9).
On the other hand, the main advantage of the simultaneous space-time finite element approach over the method of characteristics is its inherent simplicity and versatile applicability. Therefore, the present space-time finite element approach should also be well-suited for the inverse dynamics of more elaborate structures such as geometrically exact beams ¶ . Moreover, the possibility to use unstructured meshes in space-time as well as higher-order elements should be investigated in future work. The extension to shells and solids should also be possible. However, it still has to be clarified under which conditions the corresponding inverse dynamics problem is well-posed.

ACKNOWLEDGMENTS
This work was supported by the Deutsche Forschungsgemeinschaft (DFG) under Grant BE 2285/12-1. This support is gratefully acknowledged. The authors would also like to thank the two anonymous referees for their insightful comments which led to an improvement of the present work. ENDNOTES * When compared to standard (rheonomic) holonomic constraints, servo-constraints in general do not have collocation property in the sense that the corresponding constraint forces in general do not refer to the same location as the servo-constraints. † The pure Neumann problem is obtained by removing the displacement boundary conditions. ‡ As mentioned in the introduction, in contrast to servo-constraints holonomic constraints have collocation property. Concerning the imposition of essential boundary conditions in the finite element method, the constraints are commonly assumed to be holonomic. 36 § Note that the designation continuous Galerkin is attributed to continuous test functions despite the use of discontinuous weighting functions. This is in contrast to the so-called discontinuous Galerkin method which relies on test and weighting functions, both being discontinuous (cf. Reference 33). ¶ Like geometrically exact strings, the motion of geometrically exact beams is governed by second-order quasilinear hyperbolic PDEs. 22 ORCID Timo Ströhle https://orcid.org/0000-0003-0808-6315 Peter Betsch https://orcid.org/0000-0002-0596-2503