Dynamic sensitivity analysis of deployable flexible space structures with a large number of design parameters

The reliability and dynamic performance of deployable flexible space structures significantly depend on their key design parameters. Sensitivity analysis of these design parameters in the frame of multibody dynamics can serve as a powerful tool to evaluate and improve the dynamic performances of deployable space structures. Nevertheless, previous studies on the dynamic sensitivity analysis are mainly confined to planar multibody systems with a few design parameters. In this study, an efficient computational methodology is proposed to perform the dynamic sensitivity analysis of the complex deployable space structures with a large number of design parameters. Firstly, the analytical sensitivity analysis formulations of objective functions with the parameter-dependent integration bounds are deduced via the direct differentiation method and adjoint variable method. A checkpointing scheme is further introduced to assist the backward integration of the high dimensional differential algebraic equations of the adjoint variables. The flexible beams in the deployable flexible space structure are described by the locking-free three-node spatial beam elements of absolute nodal coordinate formulation. Furthermore, a parallelized automatic differentiation algorithm is proposed to efficiently evaluate the complex partial derivatives in the sensitivity analysis formulations. Finally, four numerical examples are provided to validate the accuracy and efficiency of the proposed computational methodology.


Introduction
The dynamic design of a large-scale deployable space structure is of paramount importance for a successful spacecraft mission. Due to the rapid development of both theories of multibody dynamics and computation technologies, the dynamics of such practical complex multibody systems (MBSs) can be accurately simulated [1][2][3][4].
However, the dynamic simulations merely reveal the dynamic behaviors of the present designs. These designs are usually not the optimal ones. Therefore, further explorations of the design space based on sensitivity analysis (SA) [5][6][7][8][9][10] are crucial and necessary.
The SA bridges the gap between the dynamic simulations and the gradient-based optimizations [11][12][13], optimal controls [14] and parametric identifications [15,16]. SA can be used to obtain the sensitivities of the objective functions with respect to the design parameters. The obtained objective sensitivities quantify the potential improvements of the dynamic performances of a system when certain design parameters are tweaked. Nevertheless, very few works on the SA of the flexible MBSs with a large number of design parameters are reported, since the computation efficiency is the prime bottleneck.
The simplest and most intuitive approach to perform the SA of a multibody system is the finite difference method (FDM) [17]. The implementation of the FDM includes two rounds of dynamic simulations performed by using the nominal parameters and perturbed parameters, respectively. Then, the variations of the objective functions are computed and divided by the parametric difference. Despite their simplicity, classical FDMs suffer from truncation and cancellation errors. Due to cancellation errors, an excessively small parametric step may lead to totally wrong objective sensitivities.
Thus, the SA based on analytical formulations is a promising alternative due to its superior accuracy and efficiency. The direct differentiation method [6][7][8] and the adjoint variable method [6,7,18] are the two available approaches for deriving the governing equations in an analytical form for SA.
The direct differentiation method has been widely used because of the easiness in implementation. This approach directly differentiates the equations of motion (EOM) with respect to each design parameter, producing the corresponding equations of sensitivity (EOS). The EOS can be further discretized and integrated in a forward manner along with the EOM. The generalized state sensitivities (the partial derivatives of position, velocity, acceleration and Lagrange multipliers with respect to design parameters) can be determined from the time integration of EOS. Then, the design sensitivities of the objective functions can be determined from the generalized state sensitivities by using numerical quadrature methods [17]. To minimize extra overhead [6], consistent formulations and time integrators can be used to concurrently solve the EOS and EOM due to their inherent relationship. In the computation process, however, the EOM must be differentiated with respect to each parameter separately, and each differentiation will produce a new set of EOS. Thus, this method is quite inefficient when a large number of design parameters are studied. Dopico et al. [6] proposed a direct differentiation method based on the index-3 augmented Lagrangian formulation with velocity and acceleration projections. The feature of their method is the unified handling of both holonomic and nonholonomic constraints. Bhalerao et al. [8] presented a direct differentiation method in a divide-and-conquer framework and a corresponding parallel algorithm with logarithmic complexity. However, the proposed algorithm can't be directly applied to the structures with complicated topologies, such as dense nested kinematic loops [19]. Neto et al. [9] demonstrated that the direct sensitivity analysis of the flexible multibody systems can be performed by using the mode component synthesis. Callejo et al. [20] presented a simple vehicle dynamics benchmark for sensitivity analysis based on a state-space formulation. However, the case studies showed that the above works are confined to the dynamic systems with simple topologies and a small number of design parameters.
Compared with the direct differentiation method, the adjoint variable method is more suitable for the systems with relatively few objective functions and a large number of design parameters [21,22]. The adjoint variable method enables one to evaluate the objective sensitivities directly, and avoids the calculation of generalized state sensitivities by augmenting zero terms to the objective sensitivities. The governing equations of the adjoint variable method depend on the introduced adjoint variables instead of the generalized state sensitivities. The objective sensitivities with respect to different design parameters can be computed via the common adjoint variables.
However, the governing differential equations in terms of the adjoint variables need to be integrated in backward direction. Such a feature indicates that before the backward integration starts, a large amount of extra computer memory has to dedicate to store all the state variables throughout the forward integration of the EOM. Although the checkpointing technique [23,24] can be utilized to reduce the heavy cost on computer memory, the adjoint method still consumes more memory than the direct differentiation method. 4 The adjoint variable method can be further categorized as the continuous adjoint variable method (CAVM) and the discrete adjoint variable method (DAVM) [25]. In the CAVM, the continuous EOM are differentiated prior to the temporal discretization of the time integrator. This method requires additional interpolations to obtain the state variables at the integration instants specified by backward integrators, because these instants are usually not coincident with those specified by the forward time integrator.
Yet, the DAVM proceeds in a discretize-then-differentiate manner. Thus, extra interpolations are exempted because the backward integration hits exactly the same instants specified by the forward integration. Ebrahimi et al. [12] combined the DAVM with the geometric variational integrator which preserves the symplectic-momentum during the sensitivity analysis. Schaffer [18] developed a novel piecewise adjoint method by transforming the initial value problem into the boundary value problem, which is well-suited for coarse-grained parallel computation. Pi et al. [26] performed the first-order adjoint sensitivity analysis of a flexible crank-slider described by fully parametrized beam elements of the absolute nodal coordinate formulation (ANCF).
However, according to the work by Bauchau et al. [27], the beam element used in above study suffered from locking problems. Serban et al. [28] creatively extended the adjoint sensitivity analysis to hybrid dynamic systems with memory. Sonneville et al. [29] reformulated the adjoint variable method in the frame of Lie group, which results in a condensed form of the equations. Lauß et al. [30] applied the discrete adjoint method to the identification of system parameters via the measured accelerations. Nejat et al.
[31] adopted the adjoint variable method to compute the structural sensitivity of a flexible crank-slider mechanism described by the approach of floating frame of reference. Wang et al. [32] conducted the adjoint sensitivity analysis based on the geometrically exact beam elements, and demonstrated the efficiency of the method through an optimization problem.
The afore-mentioned direct and adjoint formulations shed light on how to compute the objective sensitivities. In each time integration step of the above formulations, quite a few complex partial derivatives [18,33] are required for the design sensitivities of the mass matrix, the Jacobi and design sensitivities of the elastic forces of finite elements [34]. If they were not computed efficiently and accurately, the performance of the SA algorithm would not be successful. There are four ways [35] to compute derivatives: finite difference, analytical differentiation, symbolic differentiation and automatic differentiation (AD). The real-step finite difference method [17] used for computing 5 derivatives still suffers from the truncation and cancellation errors, thus the step size is strictly restricted and try-and-error process is required. The analytical differentiation method is accurate and efficient if the closed-form formulae can be derived and correctly implemented in the program. However, the derivation of analytical derivatives is a usually cumbersome or even impossible task sometimes, especially for the complex systems with a large number of design parameters. The symbolic differentiation can alleviate the deduction burden to some extent, but will suffer from the expression swell issues [36]. The automatic differentiation method, also called the algorithmic differentiation method, is a promising alternative compared with the above three methods. The AD method can get machine-precision derivatives efficiently with moderate efforts [37]. The method computes the derivatives by applying the chain rule of derivative calculus repeatedly to the elementary arithmetic operations in the program, without the incurring of truncation errors. The AD method can be performed in two ways, the operator overloading [38][39][40] and the source-to-source transformation [41,42], while each way also covers two modes, that is, the forward mode and the reverse mode [37]. Recent years have seen successful applications of AD to evaluating derivatives in the dynamic simulation of multibody systems [43,44], but confined to multiple rigid-body systems. In this study, the complex partial derivatives are determined from the forward mode AD based on operator-overloading.
The remaining parts of the paper are organized as follows. In Section 2, the analytical formulations of SA are deduced via the direct differentiation method and the adjoint variable method, respectively. In this study, the upper and lower integration bounds of the objective function also depend on the design parameters. In Section 3, to compute the complex derivatives of the design sensitivities and Jacobi of the nonlinear elastic forces of a beam element efficiently and accurately, a parallelized automatic differentiation algorithm is developed by using OpenMP directives. In Section 4, three numerical examples are firstly provided to validate the proposed formulations. The fourth example is to demonstrate that the proposed computation methodology can be used to the SA of deployable flexible space structures with a large number of design parameters. Finally, the main conclusions are drawn and the perspectives for future researches are outlined in Section 5.

Dynamic equations of a flexible multibody system
The dynamic equations for a flexible multibody system with m holonomic 6 constraints described by n generalized coordinates can be written as a set of differential algebraic equations (DAEs) of index-3 as follows  (2) In addition, p  a ¡ denotes the vector of p design parameters for the multibody system of concern.

SA based on the direct differentiation method
Let the dynamic response of concern be from time 0 where 1  and 2  are called the parameter-dependent kinematic criteria. In the deployable space structures, Eq. (4) can be used to describe the behaviors of the self-actuating and self-locking mechanisms [45].
Then, according to the direct differentiation method [5], the objective sensitivity can be evaluated by differentiating Eq. (3) with respect to a , as where ,

SA based on the adjoint variable method
In the above direct differentiation method, the generalized state sensitivities with respect to each design parameter have to be computed separately. This way is especially inefficient for the systems with a large number of design parameters. According to the adjoint variable method [28][29][30][31][32], the computations of generalized state sensitivities can be circumvented by a clever augment of zero terms in the objective sensitivities, and only one set of DAEs are required to be solved. 8 Using the method of integration by parts, Eq. (5) can be further cast as . From the adjoint method [5,21], the time integration of the sensitivity of the first equation in Eq. (1) can be written as where μ is the adjoint vector defined for the dynamic equations. Also, by using the method of integration by parts, Eq. (10) can be further cast as where ν is the adjoint vector defined for the constraint equations. Let Eq. (4) be differentiated with respect to vector a and then respectively multiplied by the introduced adjoint scalars 1  and 2  , the following adjoint equations can be established as Subtracting Eqs. (11), (12) and (13) from Eq. (9), the augmented equation of the objective sensitivity can be written as Then, by equating the coefficient expressions of a q and a λ in the integrands of Eq. (14) to zeros, one can obtain the following index-3 adjoint DAEs in terms of μ and ν as follows The above adjoint DAEs can be regarded as a final value problem, which can be  can be written as , .
It can be found that Eq. (16c) and Eq. (16e) fail to give a unique solution of the 10 adjoint variable 2  . To solve this problem, the following conditions should be further Again, it can be found that Eq. (16b) and Eq. (16f) are not compatible either.
Similarly, the following zero terms evaluated at moment F t are augmented to Eq. (14), where the new adjoint vectors γ and η are introduced to multiply the design sensitivities of the position-level and velocity-level constraint equations, respectively.
After adding Eq. (18) to Eq. (14), the augmented objective sensitivity can be recast as Then, the compatible equations can be established and solved for starting the backward integration of Eq. (15) as follows where 1  is obtained by solving . It is worth noting that the simplified objective sensitivities with respect to different design parameters can be computed via the shared adjoint vectors μ and ν determined from Eq. (15). Thus, the adjoint variable method is favorable for the SA of system with a large number of design parameters.

Time integration and checkpointing scheme
For the index-3 DAEs formulated by Eq. (1), the implicit Newmark integrator described in [46] is frequently used due to its simplicity and wide applicability. Taking the time integration of Eq. (7) for example, the position sensitivities where h denotes the time step,  and  are scalar parameters that affect the precision and stability of the integrator. This method is implicit and A-stable provided [47]. In the case of the backward integration required by the adjoint index-3 DAEs, the time step Applying the above integration rules to the temporal discretization of Eq. (7) at This equation is a linear one in terms of the generalized state sensitivities Eq. (15) of the adjoint variable method is solved backwards in integration time, and the backward solution process can only be performed after the forward time integration of dynamic equations is finished. To proceed the backward integration, the generalized state variables at the desired instants are required. One simple strategy is to store all the generalized state variables at the discretized time steps during the forward integration. However, this way is only acceptable for small-scale problems due to the limited memory of a computer. The checkpointing scheme proposed in [23,24] is a compromise between the storage space and the computation efficiency. Fig. 1 schematically shows the implementation of this scheme, which can be summarized as the following six steps: Computer Memories Step k Step k+1 t 0 t F Fig. 1 The backward time integration for the adjoint DAEs using the checkpointing scheme

Analytical and symbolic differentiations
The analytical differentiation is efficient and accurate, but only applicable for the systems with closed-form solutions. Computing complex derivatives often require 13 much mental work. Symbolic softwares, such as Maple and Mathematica, can be used to alleviate the deduction work of the analytical differentiation [35]. Yet, one major issue of the symbolic differentiation is the expression swell [36]. The symbolic differentiation is able to handle the control flows in the computation of complex derivatives, such as if-else statements, loops and recursions, but the produced formulae are excessively complicated. In this study, only some simple quantities such as T qq Φ λ and qq Φ q & can be determined from the analytical or symbolic differentiation. Also, it is wise to do so because these constraint-related matrices usually have regular sparse structures depending on the categories of the kinematic pairs. Once the analytical formulae of these matrices for a specific kinematic pair are derived, they can be exploited for computing the derivatives efficiently.

Real-step finite difference method
The finite-difference method (FDM) is the simplest and most intuitive technique for computing numerical derivatives. For the sake of simplicity, only the differentiation of a scalar f with respect to the parameter vector p  a ¡ is considered here.
Accordingly, the proposed formulae can be easily extended to the differentiation of vectors, matrixes and high-order tensors.
According to the Taylor series expansion, the first-order derivative of f in the real-step FDM can be approximated as where s a  is a small perturbation to the s-th design parameter of a , p s  e ¡ is the unit vector, the s-th entry of which is one and all the other entries are zeros. Despite its simplicity, the accuracy of the real-step FDM is sensitive to the perturbation steps. If the perturbation step is too large, the result is inaccurate due to truncation errors.
Otherwise, the result may be totally wrong due to canceling-off errors if the perturbation step is too small. Try-and-error is the only way to deal with this contradiction. Again, there is no universal criterion to determine the adequate perturbation step for different physical problems. Although the FDM is not used in the codes for computing the derivatives in the analytical formulations of Section 2, it is widely used for checking other methods via the adequate perturbation steps. 14

Automatic differentiation method
Compared with the above methods, the automatic differentiation method has several advantages for computing derivatives. Firstly, AD is convenient to implement because there is no need to derive the closed-form formulae of the derivatives. Secondly, AD incurs no truncation and cancelling-off errors compared with the FDM. Also, AD gets rid of the try-and-error process in the real-step FDM and enables one to obtain the derivatives with respect to the design parameters simultaneously and does not need to perturb each parameter separately.
The categories of AD can be classified into source-to-source transformation and operator overloading [37]. Source-to-source transformation works as a preprocessor to generate the codes for computing derivatives. The generated codes are usually quite complex, which prohibits an easy integration with the program of analytical SA in In the operator-overloading approach, a new data structure of AD [39] is defined to declare the independent, intermediate and dependent variables. A scalar x is extended to vector AD x as follows where x is the function value that depends on the independent variables, p x  d ¡ is the vector of the non-zero partial derivatives of x with respect to different independent variables, p is the maximal dependency number introduced to store the derivatives in a compressed format,  Indices of non-zero derivatives 3 p  Indices of non-zero derivatives P P P P P P P P P P P P P P P P P P P P P P P P 6 p    Table 1 presents the overloading rules of several common operators and intrinsic functions. As shown in the third column of the table, the three parts of the resultant are calculated separately. The first part is to perform the ordinary scalar arithmetic of the expression, the second part is to calculate the non-zero derivatives of the right-hand side based on the chain rule of differentiation, and the third part is to update the dependency indices of the right-hand side based on the union of dependency indices of the left-hand sides.
x d d n n n U Furthermore, the above-mentioned overloading process covers only the overloaded arithmetic between variables of AD data structure. The mixed mode arithmetic, like multiplication of a scalar and a variable of AD data structure, is not covered but can be accomplished similarly.
To illustrate the usage of the AD data structure and overloaded operators, consider an simple example, the AD of a scalar expression 2 . As the addition and multiplication operators are overloaded, the corresponding overloading rules in Table 1 are used such that the resulting intermediate variable in the extended form AD z can be written as x to the intermediate variable AD z . It is clearly shown that the AD data structure and overloaded operators enable the propagation of the derivatives along the computational graph automatically [37].
Moreover, the propagation will not be hindered by the control flows, such as loops, recursions, if-else statements and subroutine calls. Thus, the derivatives relay consecutively among independent variables, intermediate variables and dependent variables. The accurate derivatives of the dependent variables with respect to the independent variables can be easily extracted from the new data structure of AD.

Application of a parallelized AD to a locking-free beam element of ANCF
In this study, the flexible beam is described by the three-node spatial beam element [48] of ANCF. This element has been validated by experiment results [27]. In this subsection, a parallelized AD scheme is proposed to evaluate the design sensitivities and Jacobi of the element elastic force. The design parameters (independent variables) of the ANCF beam element with circular cross-section are assumed to be where L is the length of beam element, in R and ou R respectively denote the inner and outer radii of the cross-section. And E and  in Eq.  T  T  T  T  T  T  T  T  T  T  27  ,  ,  , ,, 22 To improve the computational efficiency, Eq. (31) can be further rewritten in a compact form as and   9 According to the work by Gerstmayr [49], the deformation gradient of a finite element of ANCF reads where 0 e denotes the vector of initial global coordinates of the nodes, and matrix  [48,50], the elasticity tensor can be firstly split into two parts, and the selective reduced integration technique can be further utilized to evaluate the element elastic forces.
As shown in Eq. (37), the evaluation of elastic forces can be considered as a volumetric integration of a polynomial function   ,, f    over the volume of the beam element. according to [51], the numerical integration can be calculated by Here, the Gauss-Legendre quadrature and Chebyshev-Gauss quadrature are used for the integrals along the beam axis and over the cross-section, respectively.

The dynamic SA of a spring-actuated five-bar mechanism
This subsection shows a benchmark for the SA of rigid body dynamics [5,6]. As shown in Fig. 6, the five-bar mechanism is actuated by two linear elastic springs without dampers, and is released from the illustrated configuration and subject to vibration under gravitational acceleration,   Fig. 7 clearly shows that the accuracy of the finite difference method strongly depends on the perturbation steps. If the perturbation step is too large or too small, the computed sensitivities will deviate from those obtained by the direct differentiation method.

The static SA of a flexible cantilever beam
To demonstrate the accuracy and efficiency of the proposed AD method in Section 3.4, the static SA of a flexible cantilever beam under a small vertical tip load is provided.
As shown in Fig. 8 The beam is meshed by the locking-free beam elements described in Section 3.4.
where the second moment of area is 3 Table 3 lists the model parameters of the two pendulums.  Prior to the sensitivity analysis of this mechanism, the performance of the lockingfree beam element is validated via the dynamic simulation of the spatial double pendulum under gravity action. Fig. 12 where F t is the termination time of concern, B X dentoes the global position of point B along X-direction. Fig. 14

The dynamic SA of a deployable seven-module antenna reflector
This subsection demonstrates that the proposed framework is able to apply to the SA of complex rigid-flexible systems. The case study of concern is a deployable sevenmodule antenna reflector for a satellite to be designed to meet the stiffness, reliability and accuracy requirements. The most important issue is that such a complicated structure has a great number of parameters in design. The initial design of such structures based on experience is far from optimal one. In this context, sensitivity analysis is quite useful for evaluating the initial design quantitatively. The case study in this subsection provides a paradigm for the sensitivity analysis of a large-scale engineering structure. O C Fig. 15 The seven-module satellite antenna reflector in the fully deployed configuration As shown in Fig. 15, the satellite antenna reflector consists of seven modules.
Module ① is located at the center and the other six modules from ② to ⑦ are surrounding module ① in a circular pattern. The deployment of each module is driven by an umbrella-like mechanism on the center rod. Moreover, each module is assembled by six identical deployable ribs radiated from the common center rod. The dimensions and topology of a basic deployable rib are also shown in the bottom of Fig. 15. The center rod, connecting rods and joint accessories are modeled as rigid bodies described by ANCF-RN [52]. The upper, lower, diagonal and vertical beams of each rib are modeled as flexible bodies meshed by the locking-free beam elements of ANCF [48].
The kinematic constraints in each module include fixed joints, revolute joints and cylindrical joints. The connections between the adjacent modules are achieved via the revolute joints. The connecting parts between the modules ③ and ④ are fixed in space, which represents a welding connection between the reflector and its satellite platform. Table 4 shows the reflector multibody model parameters. Fig. 16 shows the dynamic deployment configurations of the reflector at different simulation moments.    Fig. 17a, the obtained objective sensitivities for modules ② and ⑤ show a similar distribution, and so do modules ③ and ④, as well as modules ⑥ and ⑦. This is because these pairs of modules are symmetric about the fixed joint between the reflector and the inertial space.
It is worth noting that the finite difference method is not utilized here for comparison, because it is not practical to perturb all 672 design parameters one by one for such a complicated model. Even though the direct method and adjoint method give close objective sensitivities, their computation efficiency is quite different. As shown in Table 5, the total 33 computation time for the adjoint variable method (201163s) is about half of the one for the direct differentiation method (392727s). Thus, the adjoint variable method is advantageous over the direct differentiation method, especially for the systems with a large number of design parameters. It can also be found from Table 5 that for the adjoint variable method, although the checkpointing scheme requires the backward integration of the adjoint DAEs, the payment is necessary and worthy since the total computation time is significantly reduced compared with the direct differentiation method.

Concluding remarks
The study presents a computational framework for the sensitivity analysis of