Joint–coordinate adjoint method for optimal control of multibody systems

This paper presents a joint–coordinate adjoint method for optimal control of multi-rigid-body systems. Initially formulated as a set of differential-algebraic equations, the adjoint system is brought into a minimal form by projecting the original expressions into the joint’s motion and constraint force subspaces. Consequently, cumbersome partial derivatives corresponding to joint-space equations of motion are avoided, and the approach is algorithmically more straightforward. The analogies between the formulation of Hamilton’s equations of motion in a mixed redundant-joint set of coordinates and the necessary conditions arising from the minimization of the cost functional are demonstrated in the text. The observed parallels directly lead to the definition of a joint set of adjoint variables. Through numerical studies, the performance of the proposed approach is investigated for optimal control of a double pendulum on a cart. The results demonstrate a successful application of the joint-coordinate adjoint method. The outcome can be easily generalized to optimal control of more complex systems.


Introduction
In the design and development of many complex multibody systems, researchers have to consider trade-off between various system attributes such as sizing, performance, comfort, or cost. Computational optimization methods are almost always required for most design tasks, and the gradient information of an objective function is heavily exploited in the generation of sensitivities. To this end, a family of black box methods can be distinguished, where the gradient computation routine is agnostic of the underlying dynamics. The main approaches of this kind are finite difference method as well as far more involved automatic differentiation [1].
Conversely, one of the most accurate and computationally efficient methods of calculating derivatives of performance measure with respect to the input variables, such as design or control variables, is based on the mathematical models of multibody system (MBS). In the optimal design and optimal control (OC) of MBS, an implicit dependency exists between state and input variables. There are two major approaches to capturing this dependency: direct differentiation [2][3][4] and the adjoint methods [5,6]. Various formulations of the equations of motion (EOM) yield different adjoint systems, each characterized by different properties [7,8]. Recent works [9][10][11] have demonstrated that using constrained Hamilton's canonical equations, in which Lagrange multipliers enforce constraint equations at the velocity level, one can obtain more stable solutions for differential-algebraic equations (DAEs) as compared to acceleration-based counterparts [12,13]. This phenomenon is partially connected with the differential index reduction of the resultant Hamilton's equations.
The adjoint approach represents a comprehensive computational framework rather than a single method, where many authors develop this broad branch of research in various directions. For example, the adjoint method for solving the nonlinear optimal control problem has been applied in different engineering areas such as parameter identification [14], structural mechanics [15], time-optimal control problems [16], and feedforward control design [17]. In the conceptual setting of nonlinear mechanical systems, the control problem is investigated using the optimal control theory, where an approximate solution of the resulting set of differential-algebraic equations can be obtained employing a broad variety of numerical procedures. Furthermore, the integration routine may be combined with the discretized optimal control problem to come up with the discrete adjoint equations [18,19]. The adjoint method has also been employed in hybrid systems involving discontinuities or switching modes [20,21].
Another important aspect lies in the selection of the variables describing both the dynamic and adjoint subspaces. It is a common practice in the field of multi-rigid body dynamic computations to employ a redundant set of coordinates to describe the underlying phenomena. As a result, the dynamic equations of motion are formulated as a set of differentialalgebraic equations. Consequently, the adjoint sensitivity analysis also yields a system of DAEs that must be solved to determine the gradient of the performance measure [5][6][7][8]. On the other hand, the derivation of the adjoint system based on the joint-coordinate formulation of the underlying dynamics raises very complicated expressions for the coefficients of the resultant adjoint system. The resulting quantities are much harder to establish systematically than in the case of the redundant DAE formulation [22]. Recent works investigate the relationships and analogies between the adjoint systems formulated as DAEs and ODEs [23,24], whereas joint-coordinate adjoint formulations are also actively studied in the literature [15].
In this paper, the authors demonstrate a novel formulation of the adjoint method to efficiently compute the gradient of the cost functional that arises in optimal design or optimal control problems. Initially, a multibody system is described in terms of constrained Hamilton's equations of motion using a redundant set of coordinates. It is well known that such an approach is attractive due to its relative simplicity. However, it leads to a mixed set of differential-algebraic equations, which can be challenging to solve numerically. An alternative is to use a minimal set of coordinates equal to the number of degrees of freedom of a multibody system. Various techniques have been proposed in the literature for coordinate reduction, e.g., [25][26][27]. The same apparatus can be used to establish the relationships between redundant and independent adjoint variables. Ultimately, the adjoint system is formulated as a set of first-order ODEs with right-hand sides that can be easily presented in a closed form [28,29].
The primary importance of this research may be summarized as follows: 1. We present a novel adjoint-based method for fixed-time optimal control problems that exploits a set of independent adjoint variables incorporated in the Hamiltonian framework. 2. We show the parallels between the formulation of the equations of motion in mixed redundant-joint coordinates and the necessary conditions arising from the minimization of the cost functional. 3. We derive the adjoint system of ordinary differential equations and consistent terminal conditions expressed in a set of joint-space coordinates derived through the use of the joint's motion and constraint force subspaces. 4. We demonstrate a successful application of the proposed method to gradient computation for the optimal design of a spatial rigid body and the optimal control of a double pendulum on a cart.
This paper is organized as follows. Section 2 presents the equations of motion for rigid multibody systems formulated in redundant and joint-space sets of coordinates. Two complementary sub-spaces are introduced, which correspond to the joint's motion and constrained directions. Next, in Sect. 3.1, a design sensitivity problem is concisely recalled. The adjoint method tailored to constrained Hamilton's equations is presented in Sect. 3.2. The concept of independent adjoint variables is introduced in Sect. 4. The developed method is applied to the optimal design of a spatial pendulum and optimal control of a double pendulum on a cart. The results are shown in Sect. 5, followed by a discussion presented in Sect. 6. Ultimately, the conclusions are drawn in Sect. 7.

Hamilton's canonical equations
This section introduces the reader to the appropriate background concerning various formulations of Hamilton's canonical equations of motion. Section 2.1 is focused on the equations expressed in dependent coordinates. Subsequently, joint space formulation is derived in Sect. 2.2 by projecting the equations of motion onto the joint's motion and constraintforce subspaces, respectively. Notation and basic symbols that are used in this paper are also explained.

Constrained formulation of EOM
Let us define a set of dependent variables q ∈ R nq that uniquely describe the configuration of a multibody system. Since the number of configuration variables exceeds the number of degrees of freedom, there exist m algebraic constraints that can be symbolically written as: The time derivative of constraint equations is equal to:Φ = Φ qq , where Φ q ∈ R m×nq stands for the constraint Jacobian matrix. The time derivative of generalized coordinates,q, consists of translational (Cartesian coordinates) and rotational (e.g., Euler angles or unit quaternions) components. Conversely, this paper will assume spatial velocities v ∈ R nv as the primary variables used to formulate the equations of motion. The velocity of an arbitrary body A can be expressed as 3 and ω A ∈ R 3 refer to linear and angular velocity of a body, respectively. Moreover, spatial velocity is related to the time derivative of dependent coordinates via bidirectional, configurationdependent map:q =qB v (q)v, v = v Bq(q)q, which allows one to rewrite the derivative of the constraint equations in the following way: Here, the matrix D T = Φ qq B v ∈ R m×nv plays the role of the constraint Jacobian. Let L be the system Lagrangian function defined by L(q, v) = T (q, v) − V (q), where T = 1 2 v T Mv and V are the kinetic and potential energy of a system, respectively, and let the symbol M ∈ R nv ×nv denote a mass matrix. Canonical momenta p ∈ R nv conjugated with the velocities v are defined as: Subsequently, the Hamiltonian function of a multibody system can be expressed as H = p T v−L. Hamilton's equations of motion for constrained multi-rigid-body system can, therefore, be formulated as follows: where Q = Q ex − H T q is a sum of external non-conservative and conservative forces, respectively. The term Dλ reflects constraint reaction forces, where the quantity λ ∈ R m is a vector of Lagrange multipliers that represent constraint forces at joints distributed along the directions indicated by columns of the Jacobian matrix D T ∈ R m×nv . Now, let us augment the Lagrangian function with a term explicitly enforcing the kinematic velocity constraints (1): L * = L + σ TΦ . The quantity σ ∈ R m represents a vector of m new Lagrange multipliers associated with the velocity level constraint equations. Based on this modification, we can define the augmented momenta in the following way: It can be easily shown thatσ = λ. The physical interpretation of this relation is that σ is a vector of constraint force impulses. Moreover, from the definition of the augmented momenta follows that p * = p + Dσ. By differentiating this relation, invoking the propertẏ σ = λ, and substituting it into Eq. (3), we come up with the following formula for the derivative of augmented momenta:ṗ * = Q +Ḋσ.
Equations of motion (6a), (6b) constitute 2n v + m differential-algebraic equations of index two with position, momenta, and Lagrange multipliers as unknowns. According to the literature, this approach is often more efficient and numerically stable than the accelerationbased formulation, mainly due to the lowered differential index [10,12,13]. The formulation based on spatial velocities is preferred to the corresponding formulation presented in ref. [30] since it seamlessly integrates with the joint space formulation presented in the following subsection. This property will be even more useful when deriving the relationships between adjoint variables in Sect. 4.

Joint-space equations of motion
This section presents the derivation of joint-space Hamilton's equations via projection method. Equations of motion for open-loop system are demonstrated first and further generalized to take into account closed-loop topologies. Two orthogonal subspaces are defined next, namely joint's motion subspace H j ∈ R 6×n j dof and constrained motion subspace D j ∈ R 6×6−n j dof that can be utilized to conveniently represent joints in a multibody system [31]. For example, to describe a spherical joint interconnecting two bodies, we would define these subspaces as: where s 1c ∈ R 3 denotes a vector starting from the joint's origin to the body-fixed coordinate frame, typically a center of mass as shown in Fig. 1. The tilde operator transforms a vector to a skew-symmetric matrix, which, when multiplied by another vector, can be interpreted as a cross-product of two vectors. The relative spatial velocity of a single body can be described with the aid of the allowable motion subspace as: v rel is the joint velocity vector denoting relative velocity of body A with respect to the velocity of its parent body. Subsequently, one can show that there is a linear transformation between joint velocities β and spatial velocities v, which leads to the following global expression [9,31]: v = Hβ.
Here, both global and joint velocities are presented in a stacked vector format, whereas H(α) ∈ R nv ×k (with subindex dropped) denotes a global allowable motion subspace. The matrix H is dependent on joint coordinates α ∈ R nα , whose derivative is uniquely related to joint velocities, similarly to the absolute-coordinate case:α =αB β (α)β, β = β Bα(α)α.
Nevertheless, in many practical situations, it is correct to assume that the joint velocity is equal to the derivative of joint coordinates, i.e. β =α. Let us pinpoint that Eq. (7) is the explicit representation of the constraintsΦ = 0 [32]. Since β is a vector of unconstrained variables for open-chain systems, the following relation also holds: D T H = 0. At this point, Hamilton's equations of motion can be projected onto the global allowable motion subspace H. Premultiplying Eq. (6a) by H T while taking into account that H T Dσ = 0 leads to the first set of canonical equations: where p ∈ R k denotes joint momenta in stacked vector format, and M ∈ R k×k can be interpreted as joint-space mass matrix. Consequently, the time derivative of p can be written as: By recalling that H TḊ +Ḣ T D = d dt (H T D T ) = 0, we come up with the following set of governing equations for open-chain systems: Let us note that certain terms from global formulation, such as Q or M, also appear in Eqs. (10a), (10b). Although we have defined them as dependent on absolute coordinates q, v, it is possible to reformulate these quantities to be dependent on joint variables. Furthermore, joint coordinates can be mapped onto a vector of absolute coordinates via unique, non-linear, possibly recursive relation, that can be symbolically expressed as q = g(α).
The governing equations (10a), (10b) are applicable only to open-loop MBS. If there is one or more kinematic closed-loops in a system, the joint coordinates become dependent. A cut-joint method is usually employed to transform a closed-loop topology into an openloop mechanism [33]. Loop-closure constraint equations Θ = 0, Θ ∈ R mc are imposed to reflect that the joint coordinates are dependent. The time derivative of the loop constraints can be written as:Θ = Γ T β = 0, where Γ T ∈ R mc×n β represents the constraint Jacobian matrix. It can be shown that constrained form of equations of motion reads as: where σ c ∈ R mc is a vector of constraint impulse loads associated with loop-closure constraints.

Introduction
A prevalent task arising in the field of optimal design or control of MBS is to minimize the following performance index: where the integrand, denoted by h, represents the expression to be minimized in a fixed time horizon t f , and b ∈ R n b denotes a vector of design parameters or a set of discretized input functions. The second term in Eq. (12) is a terminal cost, suitable for prescribing a particular state of the system at the terminal time t f , which is fixed. The most efficient local optimization algorithms that are able to minimize the performance measure (12) rely heavily on the gradient of the performance measure. This computation phase is called the sensitivity analysis, and it usually constitutes a performance bottleneck. Many approaches are available to cope with the problem of efficient gradient computation. The most straightforward is a purely numerical approach that treats the underlying dynamics as a black box and computes the gradient via finite differences. However, the list of the method's drawbacks is relatively long: from high computational burden to evaluate a single derivative, through undesirable scaling in the case of multiplicity of design variables, to ambiguity in the choice of the perturbation step. Conversely, relying on the mathematical model of the dynamic equations is possible. This family of methods is significantly more difficult to implement; however, it yields the most accurate and computationally efficient results. A direct approach to the sensitivity analysis, known as direct differentiation method, employs a chain rule of differentiation to the performance measure (12): The idea is to compute the implicit state derivatives q b , v b by differentiating the underlying equation of motion (6a), (6b) with respect to the design or control variables b. This method proves to be reliable and efficient when the number of constraints is relatively large [2,3] since the state derivatives can be used for each of the constraint equations. Nevertheless, this approach requires solving dim(b) = n b sub-problems of size comparable to the underlying EOM. Concurrently, it is more often the case in the field of multibody dynamics that the number of constraints is relatively low compared with the dimension of b. The adjoint method focuses on solving different equations, thus avoiding the cumbersome computation of the state derivatives. Therefore, its efficiency is not affected by the number of design/control variables, and the total cost of computing the gradient is only proportional to the size of the underlying dynamic problem. This property makes the adjoint method the most viable option for efficient and reliable gradient computation, as will be briefly presented in Sect. 3.2.

The adjoint method
Equations of motion must be fulfilled throughout the optimization process, hence it is valid to treat them as constraints imposed on the design variables b. The key principle of the adjoint method is the inclusion of these constraints to the generic cost functional (12). For brevity, let us consider a set of DAEs describing dynamics of a multibody system, such as Eqs. (6a), (6b), in the following implicit form: F(b, y,ẏ, σ) = 0. The symbol y represents a vector of state variables of appropriate size. Consequently, the performance measure (12) may be extended to include the dynamic equation of motion in the following way: where w(t) ∈ R 2nv +m is a vector of arbitrary, time-dependent multipliers, known as adjoint, or costate, variables. Let us investigate the variation of the extended performance measure (14): The total variation δJ is a result of variations δy, δσ, δb, respectively. The variation δẏ that appears under the integral can be eliminated by integrating appropriate expression by parts: Only final time component is taken into account outside of the integral in Eq. (16) since δy| t=0 = 0 when initial conditions are explicitly prescribed. Upon substituting equation (16) into Eq. (15), we come up with the following expression: Adjoint variables have at this point arbitrary values. The goal of the adjoint method is to pick such a set of adjoint variables that would simplify equation (17) to the expression involving solely variations of design or control variables δb: Equation (18) exposes the gradient of the performance measure and allows us to obtain a unique expression for its computation by exploiting the adjoint variables: Δt -when b denotes discretized control input signals, and Δt is the discretization time-step.
To find a set of design variables or control functions that produces a stationary value of the cost functional J , the adjoint variables are required to fulfill the following set of necessary conditions: which is a general formula for the adjoint system. Equations (19a), (19b) constitute a set of DAEs that must be solved backward in time from the boundary conditions prescribed by Eqs. (19c), (19d). Equation (19d) suggests that the terminal cost S does not depend on the time derivativeẏ, which is against the assumption presented in Eq. (14). This issue can be easily circumvented; however, it involves introducing additional components that would obscure the main point presented herein. For clarity, the details on how one can treat the boundary conditions of the adjoint system are presented in Appendix A.

Adjoint in redundant set of variables
Now, let us derive the necessary conditions for the case when a rigid multibody system is described by a set of constrained Hamilton's equations. Rewriting Eq. (14) in a way that explicitly enforces the dynamic equations (6a), (6b) yields: Here, the quantities η = η(t) ∈ R nv , ξ = ξ(t) ∈ R nv , and μ = μ(t) ∈ R m constitute the entries of the adjoint variables vector w defined in Sect. 3.2. Splitting the general dynamics formula presented in Eq. (14) into more explicit equations, such as constrained Hamilton's EOM, will grant us more detailed insight into the properties of the resultant adjoint system; specifically, the relations that arise between absolute and joint-space adjoint variables. Applying the procedure described by Eqs. (14)-(19d) to the extended performance measure (20) yields the following set of equations: The coefficients A(q, v, σ) and r(q, v,v, σ, ξ, η) are quantities that raise directly from calculations presented in Sect The coefficient matrix on the LHS of equation (22b) is the same as in Eq. (6a). Equations (22a), (22b) must be integrated backward in time from known final conditions, i.e., . A more detailed discussion on the approach shown herein is presented in reference [30].

Independent adjoint variables
The goal of this section is to show the analogies that exist between constrained Hamiltonian formulation and necessary conditions for the minimization of cost functional subjected to differential-algebraic equations. As a result, a novel concept is demonstrated that allows one to introduce a set of independent adjoint variables. The key concept of this derivation lies in the relation between absolute-coordinate and joint-coordinate formulations of the equations of motion. A closer look at these relationships allows one to define a set of independent adjoint co-states. The unknown variables η, ξ introduced in equations (21a)-(21c) are dependent since Eqs. (21a)-(21c) is the adjoint system to the DAE (6a), (6b). Therefore, these quantities are referred to as absolute (i.e., constrained) adjoint variables. Conversely, the set of independent adjoint variablesη ∈ R n β andξ ∈ R n β analogous to their dependent counterparts will be referred to as joint-space adjoint variables. The joint-space dynamic equations (10a), (10b) are derived by projecting the EOM (6a), (6b) onto the global motion subspace H. Let us rewrite both sets of equations in the following way: The relation between the system of equations (23) and the set in (24) can be captured as follows: Consequently, the performance measure may be transformed in the desired way to come up with a piece of new useful information. Instead of augmenting the cost functional (12) with the absolute-coordinate equations (23), we utilize the joint-space equations (24) and invoke Table 1 Analogies between spatial velocity vector v ∈ R n v and the adjoint variable ξ ∈ R n β v ξ implicit constraint eq.
relations (25) to come up with: The variablesη(t) ∈ R n β andξ(t) ∈ R n β denote a new set of independent adjoint variables that reside in the joint-space H. By comparing equation (26) with Eq. (20), we establish the relations between joint-and absolute-coordinate co-states: The system of equations (27a), (27b) suggests thatη =ξ. In fact, this is the joint-space counterpart of the terse relation η =ξ (cf. Eq. (22a)). The simplicity of these relations comes specifically from a virtue of the Hamiltonian formulation since both momenta p * ,p and their derivatives are not premultiplied by any state-dependent coefficients. Moreover, the second equation in (23), (24) is the derivative of the first one, which is not the case, e.g., in the Lagrangian formulation. Ultimately, let us insert Eq. (21a) into Eq. (21c) and integrate the outcome by parts to get: The expression (28) can be interpreted as the implicit velocity-level constraint equations imposed on the adjoint variable ξ. Similarly, equation (27a) represents the same kinematic constraint formulated explicitly. Ultimately, equations (27a) and (28) reveal an interesting property of the adjoint variable ξ. Apparently, this quantity is somewhat analogous to the spatial velocity vector v. The relations are detailed in Table 1.
Let us now investigate algebraic constraints imposed on the physical and adjoint systems. Table 2 presents a detailed comparison of these quantities. Specifically, one can see that the adjoint constraints are shifted toward the higher-order derivative compared to the geometric constraints imposed on a multibody system. Consequently, index-1 formulation of the adjoint system (22a), (22b) involves jerk-level constraints instead of velocity-level constraints, which is the case for EOM (6a), (6b). Equations appearing in Table 2 represent implicit characterization of the constraints. An equivalent form can be obtained with the aid of the motion subspace H whose components specify the allowable velocity range space in R 6 . For example, a set of equations (27a), (27b) describes explicit velocity-and acceleration-level constraints imposed on the absolute adjoint variables. Subsequently, the differentiation of Table 2 Analogies between different algebraic constraints imposed on multibody (Φ = 0) and the adjoint system (Ψ = 0) multibody system adjoint system position-level constr.
Eq. (27b) yields the explicit form of constraints that tie up redundant and independent sets of adjoint variables:η = Hη + 2Ḣη +Ḧξ.
Once we plug Eq. (29) into (22b) and premultiply both sides by the term H T , it is possible to come up with the joint-space adjoint equations:ξ =η, A system of equations (30a), (30b) may be integrated backward in time for the unknown adjoint multipliersη andξ. The initial values can be calculated by finding the values of absolute-coordinate adjoint variables (cf. Eqs. (A.4a), (A.4c)) and projecting them onto the subspace H: Please note that the matrix H T H is invertible in non-singular configurations. Let us also pinpoint that the RHS of Eq. (30b) depends on absolute adjoint coordinates, i.e. r J = r J (q, v,v, σ, λ, η, ξ); however, the components that appear underneath are straightforward and relatively easy to formulate algorithmically [28]. The expressions that depend on the variables q, v,v, σ, λ constitute a set of parameters of the quantity r J and Eq. (30b) since their values have already been established in the forward sweep. Moreover, a system of equations (30a), (30b) constitutes a set of ODEs since the constraints are enforced explicitly via Eq. (29), and no additional algebraic constraints must be taken into account (at least for open-loop multibody systems). Once the adjoint variablesη,ξ are evaluated for the entire time domain, their absolute-coordinate counterparts can be calculated via Eqs. (27a), (27b). Ultimately, they can be utilized to efficiently and reliably compute the gradient of a performance measure (12). Let us summarize the contents of this section. The forward solution of the equations of motion yields state variables and reaction loads recorded at discrete time intervals of the temporal domain. These quantities are needed to compute the coefficients (i.e., appropriate Jacobi matrices) necessary for the adjoint sensitivity analysis. The results of almost all calculations performed in this process are to be reused for the adjoint gradient calculation (or calculated on-the-fly). Subsequently, by invoking the necessary conditions for the extremum of the performance measure (12), we came up with a set of DAEs, supplied with the relevant boundary conditions, which must be solved backward in time from the time instant t = t f to the time instant t = 0. Simultaneously, we define the dependencies between absolute and joint-space Hamilton's equations (25). These relationships allow us to derive the relations between joint-space and absolute adjoint variables (27a), (27b). Ultimately, by combining expressions (27a), (27b) with the projection methods, it was possible to develop a set of ordinary differential equations for the adjoint system formulated in terms of a set of independent adjoint variables, i.e., equations (30a), (30b). Ultimately, solving Eqs. (30a), (30b) for the adjoint multipliers allows one for the evaluation of the cost functional gradient (12). A possible approach for computation of the gradient at each step of the optimization procedure is summarized in Algorithm 1.  [30] in Eqs. (17), (27), (31).

Algorithm 1 Adjoint-based algorithm for gradient evaluation
• Define auxiliary time variable τ = t f − t , τ ∈ [0, t f ] and take into account that: d dt (·) = d dτ dτ dt (·) = − d dτ (·) • If variable-step integrator is used, exploit interpolation methods to make the grids in the forward and backward problem in agreement. 5: Recalculateξ into ξ via map in (27a) 6: Compute the gradient using an appropriate expression shown in Eq. (18) Output: gradient ∇ b J of the performance measure (12)

Numerical examples
This section presents the application results that demonstrate the validity of the methodology as well as underlying challenges encountered. We investigate the problem of open-loop control signal synthesis for a double pendulum on a cart presented in Fig. 2. The gravity acts along the negative y axis of the global frame, and a viscous friction loads are introduced at each joint in the system. The mechanism is controlled by a horizontal force applied to the cart. The vector of design parameters b consists of discretized input signals u(t), where t = 0.005 s is the discretization step. Consequently, b = [b 0 , b 1 , . . . , b N    a spline interpolation is performed to approximate it. The parameters for the system, such as masses, moments of inertia, lengths, etc., are gathered in Table 3. Although the analyzed multibody system moves in plane, it has been modeled with a spatial approach following the methodology presented in the previous sections.

Swing-up maneuver
Initial conditions for joint angles are specified to be α 1 = − π 2 , α 2 = 0 and are translated into the system stable equilibrium. The goal of this task is to stabilize the open-loop chain in the upright configuration. Mathematically, this can be described in the following way: where γ is a binary value enforcing that pendula will stop moving at the final time if γ = 1.
The initial guess is simply no actuation (u 0 (t) = 0) exerted on the cart. The gradient is calculated by solving Eqs. (30a), (30b) and conveyed to the optimization algorithm. The employed procedure utilizes ode45 Matlab function to integrate forward dynamics and the adjoint system, and fminunc to compute the desired control input. Figure 3 presents the optimization convergence rates, specifically, how the cost and first-order optimality measure progressed with respect to iteration count in the case when γ = 0. The first 30 iterations were performed with the steepest descent method, which uniformly reduced the cost to a nearzero value. Subsequently, we employed the quasi-Newton approach with the BFGS method for updating the Hessian for better convergence. In the following step, the optimization has been re-run with the last result used as the initial guess for the case γ = 1 in Eq. (32). The progress of this computation is presented in Fig. 4, while the results of both processes are displayed in Fig. 5. The plotted lines represent the absolute angle (ϕ 1 = α 1 , ϕ 2 = α 1 + α 2 ) of both pendula at the final iteration for two cases: γ = {1, 0}. Let us point out that the solid lines enter the final value with a slope, which is zero, suggesting that the final velocity is equal to zero. In summary, overall procedure consists of the following steps: 1. Set γ = 0 and use fminunc with the steepest descent method to significantly reduce the value of the performance measure (first 30 iterations in Fig. 3) 2. While γ = 0, use fminunc with the BFGS method for updating the Hessian and converge to more precise solution (iterations 31 and the following in Fig. 3; the result is denoted by dashed lines in Fig. 5). 3. Set γ = 1 and use fminunc with the BFGS method with initial guess from the former optimization (Fig. 4). The solution is recorded in Fig. 5 as solid lines.
The gradient of the objective function computed at the first iteration from ODE (30a), (30b) has been compared with its DAE counterpart (cf. Eqs. (22a), (22b)). The absolute adjoint variables can be recreated from the joint-space using a set of equations (27a), (27b).  6 The results of different gradient evaluation methods at the initial step (u(t) = 0) Figure 6 shows that these quantities have almost identical values. Moreover, we append a plot with the results taken from the complex step method [35], which additionally shows a good agreement. Lastly, Fig. 7 displays constraint violation errors at different levels: velocity-level Ψ (28), its first derivativeΨ (21c), as well asΨ. The plots refer to the ODE (cf. Eqs. (30a), (30b)) and DAE (Eqs. (22a), (22b)) formulations, respectively. The latter plot reveals that the constraint equationsΨ = 0 is fulfilled up to the machine accuracy since it is enforced explicitly in Eq. (22b). On the other hand, lower-level constraints are characterized by larger constraint violation errors, when compared against the ODE formulation. The DAE formulation exhibits the drift phenomenon due to a lack of numerical stabilization. Since the adjoint equations are solved backward in time, the most accurate results of the DAE formulation can be found at t f = 3 s, whereas the errors are accumulated as the numerical procedure approaches the initial time.

Stabilization of the cart
This test case investigates the calculation of the input control signal that will stabilize the cart while the pendula are free-falling under the gravitational force. The initial configuration of the MBS can be described in the following way: α 0 = 0, α 1 = π 4 , α 2 = 0, while the bodies' velocities are equal to zero. The optimal control approach to solving this problem requires specifying a performance measure, which can be defined in the following way:  where x := α 0 denotes an x-coordinate of the cart. On the other hand, it is possible to abandon the OC approach by adding a driving constraint equation of the form x = 0, solving only the forward dynamics problem, and recording the Lagrange multiplier λ x=0 associated with the appropriate reaction force. Let us note that the latter approach can be applied to a narrow class of problems and, e.g., it would not work in the case of Sect. 5.1. Nevertheless, the specified approach yields a very good approximation of the solution and can be utilized to compare and verify the outcome. Figure 8 depicts the input control forces calculated using the approach proposed in this paper and taking into account the driving constraint force. Consequently, Fig. 9 presents how the position of the cart changes in time for different actuation signals. The response of the cart for the optimized input signal diverges ±7 mm from the equilibrium, which is a very good fit. The starting guess for the optimization is again set to zero, i.e. u 0 (t) = 0. Figure 10 shows the convergence rate of the fminunc Matlab function, which executes quasi-Newton algorithm with a user-supplied gradient calculated from equations (30a), (30b). Figure 11 displays the gradients computed at the initial guess of the optimization and evaluated for the solution found from utilizing a driving constraint. One can see a high agreement between the methods employed. Moreover, the value of the gradient evaluated at the solution is expected to be equal to zero. This is not precisely the case in Fig. 11; however, the magnitude of the computed gradient is significantly lower when compared with the gradient evaluated at the initial guess.

Discussion
The scope of this paper is relatively broad, touching different formulations of the EOM, projection methods, and adjoint sensitivity analysis. Its main contribution lies in the novel formulation of the adjoint method based on independent adjoint variables while avoiding computing complex state derivatives, which is a necessity when standard formulation of the adjoint method is combined with the joint-space equations of motion.
In this paper, absolute-(6a), (6b) and relative-coordinate (10a), (10b) formulations of Hamilton's equations of motion are interchangeably used. By deriving the joint-space approach directly from its absolute-coordinate counterpart, the underlying dependencies between both formulations are grasped and succinctly captured by Eq. (25). These observations are further utilized to identify a set of independent adjoint variables shown in Eq. (26). The interrelations between joint-space and redundant sets of co-states are explicitly discovered in Eqs. (27a), (27b). Specifically, equation (27a) reveals a rather surprising analogy between the joint velocity β and the adjoint variableξ. Ultimately, it has been possible to formulate the adjoint system as a set of ODEs represented by a minimal set of variables by combining the derived relations with the projection methods and applying them to the absolute-coordinate adjoint system (22a), (22b). The resultant adjoint system (30a), (30b) can be solved backward in time with the aid of a standard ODE integration routine to efficiently evaluate the gradient of the performance measure.
Formulation (30a), (30b) possesses many advantages over the system (22a), (22b), which are directly inherited from the joint-space formulation (10a), (10b). The most obvious is that ODEs are generally easier to solve than DAEs and give more stable solutions during the numerical integration (cf. Fig. 7). This virtue is valid in the case of open-loop kinematic chains. The analysis of closed-loop systems requires additional efforts to include loop-closure constraints. Nevertheless, the number of algebraic constraints in joint-space framework would be significantly lower compared to the redundant setting. The scope of this paper is limited to tree-like topologies. Consequently, the number of variables that must be integrated numerically is proportional to the number of degrees of freedom, rather than the total number of redundant coordinates, which, in many practical applications, may provide a significant improvement over the absolute-coordinate formulation.
Furthermore, the coefficients and the right-hand side (r A , r B , r J ) of the adjoint system involve state-dependent jacobians of the EOM with respect to state variables. These matrices have a regular and relatively simple structure in the case of Eq. (22b), where a vast majority of expressions are readily available for implementation in closed form [28,29]. Additional details can be also found in ref. [30] (specifically in equations (17), (27), and (31)). Equation (30b) benefits from this property since its RHS is a projection of the absolute coordinate adjoint system (22a), (22b) onto the motion subspace H. Figure 12 presents various pathways one can take to generate the adjoint system. For instance, path C is advertised in ref. [30], whereas path C-D is proposed in this paper. Conversely, one can derive the adjoint system by applying variational principles to the original joint-space EOM (10a), (10b) (path A-B in Fig. 12). It can be demonstrated that such an approach produces the same set of equations as shown in Eqs. (30a), (30b), however, with the other right-hand side, i.e. r J = r J (α,α,α) [36]. Although the formulation does not involve additional Lagrange multipliers (e.g., σ, λ), the partial derivatives required to evaluate the RHS become significantly more complex and cumbersome to compute when compared with the RHS of Eq. (30b). For example, the derivative (Mβ) α (where M = H T MH) is much more involved than its absolute-coordinate counterpart (Mv) q , especially when one aims at implementing a systematic procedure.
The proposed approach is amenable to extensions that consider closed kinematic chains. The idea is based on the approach presented in ref. [33], where a subset of joint variables β * ∈ R dof is picked out of a dependent set β. The relation between β and β * is similar to that expressed in Eq. (7), i.e. β = Eβ * , which is a conversion of the velocity equation Γβ = 0. The term E ∈ R n β ×dof is the equivalent of the subspace H that is a null-space of β * . Therefore, by projecting equations (11a), (11b) onto the subspace E, we come up with the system equivalent to the one described by Eqs. (10a), (10b). Subsequently, the These areas are of current endeavors for the authors. As a final comment, the discovery of the relations demonstrated in (27a), (27b), combined with the analogies between state and adjoint variables gathered in Table 1, may pave the way to set up a fully recursive formulation for solving the adjoint system. Due to the fact that the leading matrix in the forward (Eq. (6a)) and adjoint (Eq. (22b)) problem is the same, it is expected that the Hamiltonian-based divide-and-conquer parallel algorithm [9] can be applied here to make the design sensitivity analysis performant, especially for multi-degreeof-freedom systems.

Summary and conclusions
A novel adjoint-based method that exploits a set of independent co-state variables has been developed, verified, and compared to the existing state-of-the-art formulations (redundant coordinate adjoint method, complex-step method). The joint-coordinate adjoint method is derived as an extension of the efforts presented in [30]. Parallels that exist between the formulation of Hamilton's equations of motion in mixed redundant-joint set of coordinates and the necessary conditions arising from the minimization of the cost functional are introduced in the text. The proposed unified treatment allows one for reformulation of the adjoint DAEs into a set of first-order ordinary differential equations by reusing a joint-space inertia matrix calculated in the forward dynamics problem. The approach is relatively simple to implement since the core partial derivatives are calculated with respect to the redundant set of coordinates and projected back onto the appropriate joint's motion and constraint-force subspaces. The validity and properties of the approach are demonstrated based on the optimal control of a double pendulum on a cart. Through numerical studies, the performance of the proposed approach has been quantified, especially in terms of constraint violation errors resulting from backward integration of the adjoint system. Applications show that the joint-coordinate adjoint method is stable, and the accuracy of computations is higher as compared against redundant-coordinate counterpart. Although the method has been tested only in two scenarios, the proposed approach is applicable to optimal control of multibody systems larger than those presented in the text within longer time horizons.

Appendix A: Boundary conditions for the adjoint variables
This appendix shows the derivation of terminal conditions that arise in the optimal control problem presented in this paper. In the first part, the formulation is general. Subsequently, it is tailored to the Hamiltonian approach. Section 3.2 embraces a formulation of a general method for the adjoint system. The necessary conditions for the adjoint variables include the DAEs (19a)-(19b) as well as the boundary conditions for the adjoint multipliers. Equations (19a)-(19d) allows only to prescribe the terminal configuration of the system, which is dictated by the fact that Sẏ must be equal to zero to be compatible with Eq. (19d). This implies that there is no dependency between terminal cost and the termẏ, which is clearly against the assumption specified earlier in Eq. (14). This issue can be resolved by defining additional adjoint variable x ∈ R 2nv +m at the terminal time t f and appending the quantity x T F| t f = 0 to equation (14). Its variation reads as: which transforms equations (19c)-(19d) into the following form: The term F T b x| t f must be added to the resultant variation of the performance measure (18), hence the formulas for the gradient should be updated accordingly.
In the case of Hamilton's formulation, the LHS of Eq. (A.1) expands to the following expression: where ν ∈ R nv and ε ∈ R m is a subset of the more general variable x. By computing the variation (A.3) and appending it to Eq. (20), we obtain the following boundary conditions for the adjoint variables: The term r where the compact and systematic structure is achieved by defining some additional quantities: (B.9) One can associate the tilde symbol in the Eqs. (B.9) as an indication of the explicit partial derivative of the equations of motion (6a), (6b) with respect to the appropriate variable, implied by the lower index. Equations (21a)-(21c) constitute a set of necessary conditions for expressing the variation of the cost functional solely in terms of the variation of the design variables. Upon solving these equations and determining the unique adjoint variables, equation (B.8) significantly simplifies into the analog of Eq. (18): where ṗ * b = Q b based on the earlier assumption. Equation (B.10) can be used to directly compute the gradient of the performance measure.