A frequency domain approach for parameter identification in multibody dynamics

The adjoint method shows an efficient way to incorporate inverse dynamics to engineering multibody applications, as, e.g., parameter identification. In case of the identification of parameters in oscillating multibody systems, a combination of Fourier analysis and the adjoint method is an obvious and promising approach. The present paper shows the adjoint method including adjoint Fourier coefficients for the parameter identification of the amplitude response of oscillations. Two examples show the potential and efficiency of the proposed method in multibody dynamics.


Introduction
Applications of the adjoint method to solve a variety of optimization problems in engineering sciences are widespread. Much attention to this approach has been paid recently in the context of multibody systems (see, e.g., [1][2][3][4][5][6][7]) in the field of optimal control, sensitivity analysis, and parameter identification. In [8], the adjoint method is seen as a special case of linear duality, which dramatically improves the efficiency of the computation only solving the dual problem. The basic idea of the adjoint method, e.g., as presented in Nachbagauer et al. [7], is the enhancement of the cost function by the system equations of motion including specific system parameters or controls to identify. By including the system equations of motion in the cost function, adjoint variables have to be introduced, leading to the dual problem when solving for these adjoint variables. Minimizing the cost function leads to a B S. Oberpeilsteiner stefan.oberpeilsteiner@fh-wels.at 1 classical optimization problem identifying unknown parameters of the system, as, e.g., the mass or inertia parameters of a body or stiffness and damping parameters of a spring-damper force involved. These prescribed examples lead to an identification of parameters in time domain. Various applications of the adjoint method in multibody dynamics for optimal control problems and parameter identifications in time domain can be found in recent works, e.g., in [7]. Adjoint sensitivities have been used in a penalty formulation in time domain for a full vehicle model in [9].
Moreover, from the experiment point of view, the frequency domain plays a key role when analyzing complex multibody systems. Very often only frequency ranges can be investigated in detail. Too low frequencies may not be measured, e.g., with acceleration sensors, and too high frequencies are mainly caused by measuring noise. Identification in the time domain would lead to some kind of best-fit solution. Hence, the goal of the identification is in general to fit a special frequency range. In [10], a system identification for vehicle dynamic applications has been presented based on impulse-momentum equations using a transfer function written as a frequency response function in order to take into account low and high frequency ranges. Spectral element techniques for parameter identification can also be found in the field of layered media in structural dynamics [11]. Therein, the characteristic function of the system, combining the response and impulse force function of the system, is represented in the frequency domain. The transfer function which characterizes the system in the frequency domain is then given as a Fourier transformation. The wavelet transform is used in [12] as a time-frequency representation for the determination of modal parameters of a vibrating system. Therein, natural frequencies, damping ratios, and mode shapes are estimated in the time domain from output data only. A wavelet-based approach for parameter identification is as well presented in [13]. Systems with cubic nonlinearities and systems undergoing both continuous and stick/slip motion have been addressed therein.
The latter mentioned works emphasize the importance of the spectral analysis of the system in order to understand the behavior of the system and consequently be capable of efficient parameter identification. The present paper shows a method for parameter identification in complex multibody systems in the frequency domain. A combination of the adjoint method and classical Fourier analysis for the identification of the amplitude response is presented herein as a novel approach and is applied to engineering problems.

Problem definition: cost function in terms of Fourier coefficients
Let us consider the system equations of motion in first order forṁ where p may describe the unknown parameters of the system. For the sake of simplicity, we assume the system has only one output which depends on the states y(t) := g(x). By applying classical Fourier analysis, y(t) can be approximated by in which ω k represents the angular frequency of the kth harmonics, each of which is assigned to the appropriate value of its amplitude A 2 k + B 2 k . The corresponding Fourier coefficients A k and B k are defined by For example, it might be of interest to determine a set of parameters p in such a way that the measured amplitudes Ā 2 k +B 2 k , k = 1, . . . , N of the kth harmonics are best approximated by the amplitudes of the system. Therefore, the goal is to find p such that an error function of the form is minimized. However, the problem may as well be specified in the so-called Mayer form [14]. For that purpose, the Fourier coefficients are introduced by the differential equationsȧ with the corresponding initial conditions a k (0) = b k (0) = 0, yielding A k = a k (T ) and B k = b k (T ). Hence the cost function is considered as a function of the final values of a k and b k , i.e., J = J (A k , B k ).

The adjoint gradient computation
Following the basic idea presented in Nachbagauer et al. [7], the adjoint method is applied to the cost function in Eq. (5). In a first step, the cost function is enhanced by the system equations in Eqs. (1), (6), and (7), leading to Herein ξ represents the vector of adjoint variables corresponding to the state vector. Moreover, α k and β k with k = 1, . . . , N are the adjoints related to the Fourier coefficients. Note that ξ (t), α k (t), and β k (t) can be arbitrary time functions at this point, since J = J , if Eqs. (1), (6), and (7) are satisfied. Let us now consider a forward solution x(t) of the system equations in Eq. (1) for a set of parameters p. A variation δp will result in the variations δx(t), δa k and δb k , respectively, and moreover, in a variation δJ . Up to a first order approximation, δJ is given by After applying integration by parts to the terms including δẋ, δȧ k (t), and δḃ k (t), the variation can be written in the form where the initial conditions δx(0) = 0, δa k (0) = δb k (0) = 0 and end conditions δa k (T ) = δA k , δb k (T ) = δB k are used. The computation of δx, δa k , and δb k can be circumvented, if the factors multiplied vanish. First, the terms including δa k and δb k disappear, iḟ α k =β k = 0, i.e., if α k = const. = α k (T ) and β k = const. = β k (T ). Second, the terms including δx vanish, if the adjoints ξ are defined bẏ The boundary conditions are chosen such that ξ (T ) = 0 in order to eliminate the coefficients of δx(T ) in Eq. (10). Finally, the terms multiplied with δA k and δB k can be eliminated by defining α k (T ) and β k (T ) from With x(t) from the forward solution of the system equations in Eq. (1) and the backward solution for ξ (t) of the adjoint system in Eq. (11), the variation of J is readily given by which is in accordance with the variation of the original cost function δJ . Hence, the gradient of J reads

Application to multibody systems
In most cases, multibody systems are described by a system of differential algebraic equations (DAE) in which q denotes the vector of redundant generalized coordinates, M the symmetric mass matrix, and f the vector of generalized and gyroscopic forces. Due to the algebraic constraints C(q) = 0, the equations of motion are extended by constraint forces of the form C T q λ, where C q represents the constraint Jacobian and the vector λ includes the Lagrange multipliers. Moreover, p is the vector of system parameters.
Using the additional variables v =q, the equations of motion can be reformulated as the following first order system of equations:q In this setting, the system output depending on q and v is denoted again by y(t) = g(q, v). The differential equations for the Fourier coefficients can be formulated analogously to Eqs. (6) and (7) byȧ Again, the goal is to find a set of parameters p such that the cost function in Eq. (5) is minimized. The enhancement of the cost function by the system equations in first order form in Eq. (17) and by the additional differential equations for the Fourier coefficients in

Eqs. (18) and (19) leads to
in which ξ and ζ represent the adjoint variables corresponding to the multibody system, μ pertains to the constraint equation, and α k and β k with k = 1, . . . , N are the adjoints corresponding to the Fourier coefficients a k and b k , respectively. At this point, the variables ξ (t), ζ (t), μ(t), α k and β k may be chosen arbitrarily. The variation of the function J is given by Integrating by parts of the terms with δq, δv, δȧ k , δḃ k and assuming that δq(0) = 0, δv(0) = 0, δa k (0) = 0 and δb k (0) = 0, as a consequence of prescribed initial conditions for q, v, a k and b k , yields The computation of the variations δa k , δb k , δq, δv, and δλ can be circumvented, if the factors multiplied vanish. In case of the adjoints α k and β k , constant values α k = α k (T ) = ∂J ∂A k and β k = β k (T ) = ∂J ∂B k are used to fulfillα k =β k = 0. The adjoint variables ξ (t), ζ (t), and μ(t) have to be chosen such that the adjoint equationṡ are used. At this point the adjoints can be chosen arbitrarily at t = T . For the sake of simplicity, they are set to zero, ξ (T ) = ζ (T ) = μ(T ) = 0, in order to eliminate the corresponding boundary terms.
It has to be mentioned here that the symmetry of the mass matrix M = M T has been used. If Eqs. (23) are satisfied, the variation δJ reduces to which directly relates the independent variation δp to the variation of the objective function. Hence, we may conclude that the gradient of J is given by 5 Numerical examples

Cart pendulum system
In order to present the performance of the identification method, we study a system of pendula connected to a cart performing a one-dimensional motion. As it can be seen in Fig. 1(a), the pendula are interconnected with rotational springs, and therefore this configuration represents a discretization of a rotating flexible beam. In this example we assume the parameters of the flexible beam, the stiffness c f , and damping coefficient d f to be unknown. The parameter d c , representing the cart's friction, remains untouched during the identification process at a prescribed value. The desired spectrum needed for the computation of the cost function J in Eq. (5) is generated by using the angle ϕ 1 (t) as a system output. For this purpose a numerical simulation utilizing the parameters listed in the table in Fig. 1(b) is performed in order to obtain some kind of virtual measurement. The resulting amplitude spectrum is shown in Fig. 2.
The actual parameter identification is initiated with the parameters c f = 8.5 Nm and d f = 0.15 Nms, and therefore the initial spectrum presented in Fig. 3 differs from the desired one. In order to include an additional model uncertainty, the prescribed damping parameter d c is modified to 0.1 Ns/m for the identification of c f and d f , whereas the value d c = 0.01 Ns/m is used for generating the virtual measurement. This allows for pointing out the main advantage of the presented method, which is the possibility to filter data and  In order to demonstrate the difference from identification performed in the time domain, in Fig. 4 the pendulum angle ϕ 1 is plotted for both parameter combinations. Due to the increased damping coefficient d c , the amplitude of the initial trajectory differs significantly from the measured one, especially at the end of the experiment. In the time domain there is no opportunity of filtering data, and therefore this approach would lead to incorrect parameters c f and d f .
By using a quasi-Newton method like the BFGS algorithm for finding a minimum of J and incorporating the computed gradient of Eq. (25), the solution can be found within 10 iterations. In Fig. 5(a) the convergence history for the optimization process is shown. The contour plot in Fig. 5(b) gives an impression of the complex shape of J (c f , d f ) and displays the optimization path taken by the BFGS algorithm. The final parameters gained by using

Identification of torsional vibration damper (TVD) parameters
In this section the presented method is applied to a model of a four-cylinder engine, shown in Fig. 6. The detailed model equations are given in Appendix A. The goal is to identify the parameters of the engine's torsional vibration damper (TVD), which is described by two Maxwell elements. The TVD is installed in order to reduce torsional oscillations of the crankshaft, which show large amplitudes at the 6th engine order.

Model structure
Crankshaft The torsional vibration modes of the crankshaft are crucial for the parameter identification process. Therefore, six lumped masses resulting in six degrees of freedom q 1 , . . . , q 6 (see Fig. 7  Piston In the used model, each piston features only one translational degree of freedom. The mass of each piston is denoted as m p and the piston's effective area as A P .

Dual mass flywheel
The primary side of the dual mass flywheel (DMF) is mounted on the right end of the crankshaft (see Fig. 7). Hence, its moment of inertia is assigned to q 1 . Instead of introducing a degree of freedom for the secondary side of the DMF, the prescribed angle q runup (t) is used. A nonlinear torsional spring and a linear damping element is used for connecting the primary with the secondary side.

Torsional vibration damper
The torsional vibration damper (TVD) is installed for reducing the internal torsional vibrations of the crankshaft. Within the case of the TVD, a flywheel ring is gliding in a viscous fluid. Usually, the mathematical model of the TVD is approximated by Maxwell elements. In this contribution, two Maxwell elements and one parallel spring are used. In Fig. 9 a schematic description of the model and the degrees of freedom used are presented. The accuracy of the parameters c * 1 , c * 2 , d * 1 , d * 2 , and c * par supplied by the manufacturer is not satisfactory, and therefore the values of the four parameters c * 1 , c * 2 , d * 1 , and d * 2 are adjusted by using parameter identification in the frequency domain. The masses m * 1 and m * 2 are set to zero. Pulley wheel The pulley wheel used for driving additional aggregates introduces another degree of freedom (q 26 ), which is connected to the TVD hub using a linear spring/damper with parameters c P W and d P W .

Cylinder pressures
The cylinder pressure is given by a two-dimensional map depending on the rotational speed and the crankshaft angle. The pressure is applied on each piston in accordance with the firing order.
Run-up of the engine In order to simulate the run-up performed on the real test bench, the rotational speed of the secondary side of the DMF is increased up to the final rotational speed. The ramp used for the simulation is given in Fig. 10. The run-up to the final rotational speed itself is done within the time interval [t 0 , t 1 ].

Results of the parameter identification
In this section the identification of four parameters c * 1 , c * 2 , d * 1 , and d * 2 of the torsional vibration damper is presented. The main purpose of the TVD is to reduce torsional vibrations of the flexible crankshaft caused by the periodic and dynamic loads (e.g., cylinder pressures). Hence, the twist angle of the crankshaft y(t) = q 6 (t) − q 1 (t) is chosen to be transformed into the frequency domain. As the measurements performed on a test bench commonly result in spectra for different engine orders, they are used in this investigation, too. Basically, an engine order relates the Fourier coefficients with the rotating frequencies of the engine's crankshaft. In case of the four cylinder engine, the amplitude of the 6th order is dominated by the parameter of the torsional vibration damper. In contrast to Eqs. (6) and (7), the differential equations for computing the Fourier coefficients arė where η(t) represents a window function and ω k the kth frequency of interest. It has been shown that the Hanning window (or Von-Hann window) given by is a good choice for filtering the system output. The upper time limit t u and the lower time limit t l determine the borders of the window function. Note that an amplitude correction factor g c = 2 is required (for further details see Appendix B). In Fig. 11 the system output y(t) is shown for the entire time interval. Moreover, a small interval t ∈ [t l , t u ] of y(t) is depicted in detail. The black line shows the original system output y(t), while the dashed line shows the Hanning window function η(t). The blue line is the multiplication of y(t) with η(t), which is used for the Fourier transformation. Due to the slowly increasing ramp shown in Fig. 10, assuming a steady state with constant angular velocity is valid. For the rotational speed n the time interval [t l , t u ] is given by Here two periods of the base frequency, which is 2/n for a four stroke engine, are contained in the Hanning window (for further details see Appendix B). In Fig. 12 the dashed line shows the vibration angle corresponding to the 6th engine order using initial parameters of the TVD. The dotted-dashed line (green) represents the measured vibration angle of the 6th engine order on a test bench. The black line (with the triangle symbol) shows the vibration angle of the 6th engine order with the identified TVD parameters. In the interval n ∈ [n 0 ,n 1 ], highlighted in Fig. 12, the deviation between the measured amplitude and the simulated amplitude is included in the cost function of Eq. (5). Hence, a significant improvement can be seen in this range compared to the simulation utilizing the initial parameters of the torsional vibration damper.
In Fig. 13 the engine orders 2, 4, and 8 are depicted, too. On the one hand, the engine orders of the simulation with the initial parameters and, on the other hand, the engine orders of a simulation with the identified parameters of the TVD are shown. As expected, the spectra of the orders other than the 6th are only slightly affected.

Conclusions
The present paper shows a new method for parameter identification using the amplitude response of oscillations in multibody system dynamics. The proposed method combines the  classical Fourier analysis with the adjoint sensitivity analysis for the gradient computation in an optimization problem. Using the spectrum of a system output helps in understanding the behavior of a system whereas the interpretation of the time domain data is not always promising because of excessive time ranges, noisy signals, or systematic errors.
The method is applied to a nonlinear cart-pendulum system. In order to demonstrate that the method copes with system uncertainties in the example, one parameter, which is not part of the identification, is set to a value different from the one used for generating the measurement. Therefore, the desired spectrum cannot be obtained to its full extent, but only in the given frequency range. In the time domain this would lead to an optimization problem that is not well-posed.
Further, an identification is performed using the model of a four-cylinder engine. Here, the robustness of the method and also the capability to deal with larger systems are pointed out.

Appendix A: The model equations
The equations of motion of the four-cylinder engine are summarized here. As mentioned in Sect. 4, the model equations are in the form M(q)q + C T q (q)λ = f(q,q, p, t), where q is the vector of redundant degrees of freedom given by Note that we have canceled the degree of freedom of the secondary side of the dual mass flywheel, because its motion is given by a function q runup (t). Hence, we exchange the degree of freedom of the dual mass flywheel q 23 by the given function q runup (t).
The corresponding force vector is where the 6th entry is The first term in the force vector F DMF (q, q 0 (t)) is given by a nonlinear stiffness characteristic.
The crankshaft, conrods, and the pistons are described with redundant coordinates. Hence, we introduce four algebraic constraint equations for each cylinder, where r cs is the radius of the crankshaft, l cr is the length of the conrod, and s cr is the distance to the center of mass.
Finally, the vector of constraint equations C(q) is given by r cs sin(q 2 ) − l cr sin(q 9 ) q 7 + r cs sin(q 2 ) − s cr sin(q 9 ) q 8 − r cs cos(q 2 ) − s cr cos(q 9 ) q 19 − r cs cos(q 2 ) − l cr cos(q 9 ) r cs sin(q 3 ) − l cr sin(q 12 ) q 10 + r cs sin(q 3 ) − s cr sin(q 12 ) q 11 − r cs cos(q 3 ) − s cr cos(q 12 ) q 20 − r cs cos(q 3 ) − l cr cos(q 12 ) r cs sin(q 4 ) − l cr sin(q 15 ) q 13 + r cs sin(q 4 ) − s cr sin(q 15 ) q 14 − r cs cos(q 4 ) − s cr cos(q 15 ) q 21 − r cs cos(q 4 ) − l cr cos(q 15 ) r cs sin(q 5 ) − l cr sin(q 18 ) q 16 + r cs sin(q 5 ) − s cr sin(q 18 ) q 17 − r cs cos(q 5 ) − s cr cos(q 18 ) q 22 − r cs cos(q 5 ) − l cr cos(q 18 ) with ω k = kω = k 2π T and the period duration T . Using the Hanning window function the Fourier-coefficient of the ith engine order (i ≥ 1) can be computed bȳ where n is the number of periods to be considered. For n = 1 the Fourier coefficientā i is and for n > 1 the Fourier coefficientā i is It has been shown that more than one period is required when using a Hanning window. Otherwise, the amplitudeā i includes information from other harmonics, which leads to wrong results. Furthermore, the amplitudeā i has to be corrected by a factor of 2 to be comparable with measured amplitudes.