Nonlinear Damping Identification in Nonlinear Dynamic System Based on Stochastic Inverse Approach

The nonlinear model is crucial to prepare, supervise, and analyze mechanical system. In this paper, a new nonparametric and output-only identification procedure for nonlinear damping is studied. By introducing the concept of the stochastic state space, we formulate a stochastic inverse problem for a nonlinear damping. The solution of the stochastic inverse problem is designed as probabilistic expression via the hierarchical Bayesian formulation by considering various uncertainties such as the information insufficiency in parameter of interests or errors in measurement. The probability space is estimated using Markov chain Monte Carlo MCMC . The applicability of the proposed method is demonstrated through numerical experiment and particular application to a realistic problem related to ship roll motion.


Introduction
The energy dissipation in dynamic system is often neglected or overly simplified in engineering design.However, there are many practical cases where an appropriate nonlinear model of damping is essential and an effective identification is needed.One example could be a material with high damping which is widely used for vibration and noise reduction.In addition, the increasing demands for enhanced and reliable performance of vibrating structures are requiring appropriate modeling of nonlinear damping for mechanical system.This is because of the fact that different types of nonlinear model could yield different responses of a system.
The identification problems of nonlinear system are very important in engineering in order to avoid unwanted instabilities or failure of the system.There have been considerable numbers of studies on the identification of nonlinear damping.For example, Iourtchenko et al. 1 successfully proposed an identification method.Iourtchenko and Dimentberg 2 further described a procedure based on the stochastic averaging method for inservice identification of the damping characteristic from a measured stationary response.Mohammad et al. 3 introduced a method to estimate nonlinear damping for linear and nonlinear structures.Tomme 4 introduced a method, based on modal analysis, to evaluate damping from measurements in materials and structures.The related studies mentioned previously are parametric identification methods, which focus on direct estimations of coefficients of assumed form for nonlinear damping in the nonlinear system.
Recently, Jang et al. 5 proposed a nonparametric identification method based on deterministic inverse approach, Landweber regularization method.Jang et al. 6 proved the applicability of the method through experimental work.The results presented in the literatures 5, 6 are very encouraging because this method does not require any assumption of the form of nonlinear damping unlike parametric identification method.However, this method has some limitations because it is dependent on the deterministic inverse approach.This kind of technique is content with singling out one solution without quantifying the related uncertainties or rigorously considering the stochastic nature of data noise driven by approximation errors, rounding errors, and measurement errors.Furthermore, it is not possible to explain the variability of the possible solutions because there exists an ensemble of inverse solutions consistent with the data due to unstable nature.
In this study, we proposed a new method based on a stochastic inverse approach to identify a nonlinear damping of nonlinear oscillatory system.The unique features of this paper are as follows.Firstly, an original method for the identification of nonlinear damping is proposed.The method can also be classified as nonparametric method which does not require any assumption of the form of nonlinear damping.That is, this method can be applicable to a system when no or little information is available about the actual form of nonlinear damping unlike the parametric method 1-4 which is limited to be used when there is enough knowledge about the nonlinear damping.Secondly, we formulate a stochastic inverse problem of identifying nonlinear damping by introducing a stochastic state space.The stochastic inverse problem studied here can naturally resolve the computational difficulty of the deterministic inverse model 5-10 , that is, lack of stability.In addition, it can quantify various uncertainties arising from an insufficiency of information on parameter of interests or measurement errors.Thirdly, we develop methodologies enabling probabilistic modeling for the solution based on hierarchical Bayesian formulation 11-14 .There have been lots of studies on the use of the hierarchical models applied to very challenging problems which arise in atmospheric chemistry, environmental sciences, or metabolic models 15-17 .Although this study considers a relatively simple mathematical model which is expressed as nonlinear single degree of freedom system, the present work may provide some insights into the work involved in nonlinear system identification.The main novelty is that we do not focus directly on the hierarchical model, but rather, on the identification procedure.Lastly, the method presented in this paper only requires system output, that is, motion responses.If the input, the prescribed external force, is acted on the system, then the nonlinearity of the system is easily found through the "black-box" model.
The outline of this paper is as follows.Section 2 introduces the mathematical description for the problem.In Section 3, nonlinear damping identification problem is discussed with its instable characteristics.Section 4 illustrates the formulation of stochastic inverse problem.In addition, probabilistic modeling of the stochastic inverse solution with the explanation of Markov random field 18, 19 is also described.Hierarchical Bayesian formulation and exploration of the posterior state space via Markov chain Monte Carlo 20 are introduced in Sections 5 and 6, respectively.Section 7 illustrates the identification procedure with numerical examples.The particular application to a realistic problem is demonstrated in Section 8.The problem is concerned with identifying nonlinear roll damping of a ship from free-roll decay experiment.Finally, we summarize the result and make conclusion in Section 9.

Mathematical Description
where ξ 1 and ξ 2 satisfy and Wronskian W is

Nonlinear Damping Identification
The focus of this study is to identify nonlinear damping of the nonlinear dynamic system when dynamic responses are observed.The identification of this unknown nonlinear damping becomes feasible with measurement of responses during 0 < t < T. Let y denote the measured response data, that is, y y 1 , y 2 , . . ., y m ; then we can obtain the following system for unknown nonlinear damping θ τ j B y τ j , ẏ τ j : where The identification of the unknown nonlinear damping B can be achieved by solving 3.1 regardless of any assumption on the form of nonlinear damping.Therefore, this method can be classified as the nonparametric method.The merit of nonparametric method is that this method does not require the specific form of a system explicitly.
However, unfortunately, 3.1 becomes ill-conditioned system causing computational difficulty to obtain inverse solution.This is mainly because the system 3.1 is obtained by approximating the first-kind integral operator.According to inverse problem theory 22, 23 , the approximation of the first-kind integral operator yields ill-conditioned system regardless of the choice of approximate methods such as the quadrature method and the Galerkin method 22 with orthonormal basis functions.Thus, the identification of nonlinear damping is very difficult to achieve both mathematically and numerically because even very small amount of noisy data in the measured displacement can be amplified and can give an effect on an inverse solution.

Stochastic Inverse Formalism
To address the difficulty discussed in the previous section, we newly formulate a stochastic inverse problem for the identification of the nonlinear damping by introducing a stochastic state space.A stochastic inverse problem can be formulated by considering variables g and θ in 3.1 as random variables.Let us define a multivariate random variables Θ : Ω → R m ; then a multivariate random variable G : Ω → R n can be defined by G HΘ.

4.1
For each realization g {g 1 , g 2 , . . ., g n } of g G, the following relationship is given: where θ {θ 1 , θ 2 , . . ., θ m } is a realization of the random variable Θ. Equation 4.2 can be interpreted as a stochastic function corresponding to 3.1 because it is a function of random variables.Therefore, the problem of finding θ given a realization g is the stochastic inverse problem.Unlike the ill-conditioned system in 3.1 , 4.2 becomes wellposed in an expanded stochastic space and gives a solution with distribution of random unknowns.
Assume that the dynamic response y is observed, then a single realization of directly observable parameter g can be obtained from 3.3 , and then the solution to the stochastic inverse problem can be expressed by Bayes' rule 11-14 : where p g | θ , p θ , and p g are known as the likelihood, the prior probability density function, and the normalizing constant, respectively.It is noted that the normalizing constant is not necessary to be computed for sampling procedure 11-14 .Posterior probability density function p θ | g in 4.3 can then be evaluated as follows: It is usually easy to obtain the likelihood p g | θ provided that the measurement data is contaminated by Gaussian random noise ε with zero-mean and standard deviation σ.The likelihood then depends on the distribution of random noise.Consequently, the likelihood is given by where • 2 and m refer to Euclidean norm and the number of measurements, respectively.Next, we should consider the prior probability density function p θ which reflects the information of the unknown θ prior to collecting the data.In this paper, we adopt Gaussian model with Markov random field 18, 19 , which is the most popular model for prior.A Gaussian Markov random field for the unknown θ is of the following form: Here, n is dimension of θ, λ is the scaling parameter, and the matrix W ∈ R n×n is determined as where n i is the number of neighbors for the point i.Then, with the likelihood 4.5 and the prior distribution 4.6 , posterior probability density function 4.3 can be formulated as

Hierarchical Bayesian Formulation
The maximum a posteriori MAP estimate 11-14 for posterior probability density function 4.8 can be derived as and the MAP estimation is similar to Tikhonov regularization method 22-24 , the representative deterministic inverse approach, which minimizes the functional with η λσ 2 .Therefore, it can be easily found that the scaling parameter λ and the standard deviation σ play an important role like the regularization parameter in the deterministic approaches.Furthermore, the standard deviation σ describes the noise level of the measurement data, and it is hard to quantify directly, except when the experiment for acquisition of data is not carried out repeatedly.In almost all deterministic approaches, the performance of methods highly depends on the regularization parameter and this parameter needs to be chosen a priori 5, 6, 22-24 .However, the choice of the optimal regularization parameter has never been trivial.
Unlike the deterministic inverse approach, a hierarchical Bayesian model, which is known as modern Bayesian analysis, provides a flexible way to choose hyperparameters σ and λ by treating these parameters as random variables 11-14 .As a result, the hyperparameters are automatically chosen through the hierarchical Bayesian formulation and there is no need to worry about the choice of optimal regularization parameter unlike deterministic inverse approach.This is one of distinct features of stochastic inverse problem.The hierarchical Bayesian model 11-14 for the present identification problem can then be formulated as

Posterior State Exploration
The posterior probability density function 5.3 corresponding to the solution of inverse problem in 4.2 was formulated.This probability function quantifies associated uncertainties and includes a stochastic nature of data driven by various errors.The posterior probability density function 5.3 needs to be estimated for the purpose of identification of nonlinear damping which is desired solution of the stochastic inverse problem.The estimation can be achieved through Markov chain Monte Carlo combined with Metropolis-Hastings and Gibbs algorithms 25, 26 .
In Algorithm 1, N MCMC is the total number of samples, is an easy-to-sample proposal distribution, and the conditional p θ i | θ −i , λ, σ which follows a multivariate Gaussian distribution can be expressed by In a similar manner, the full conditional p λ | θ, σ and p σ | θ, λ can also be derived as 6.2

Numerical Experiments
In order to demonstrate the workability of the method proposed to identify nonlinear damping, the numerical experiments are carried out.As a numerical example, Van der Pol Equation 27 is considered as where μ 1.Comparison of 2.1 and 7.1 shows that B y, ẏ μ 1 − y 2 ẏ and M y y.As first step toward the identification of nonlinear damping, the nonlinear model 7.1 is simulated using numerical integration scheme such as the Runge-Kutta method with the initial conditions y 0 1 and ẏ 0 0. Figure 1 plots the simulated responses and phase diagram with respect to time.In this study, the simulated displacement is considered as the measured data for the inverse identification of nonlinear damping.If the displacement is measured, the pseudo-displacement g which is necessary data for the identification of the nonlinear damping can be calculated through 3.3 and is plotted in Figure 2. In practice, the measured data is contaminated with the various noises such as approximation errors and rounding errors.In order to examine this kind of uncertainties, the synthetic noisy data is generated by adding Gaussian random noise with zero mean and unit standard deviation.The noisy data g δ g e with 10.1% noise is also plotted in Figure 2.
To show the instability of the original mathematical model in 3.1 , the numerical solution without any numerical treatments is first shown in Figure 3.That is, after decomposing 3.1 into singular systems 28 with singular values and vectors, the pseudoinverse 22, 23 is conducted.For the numerical experiments, the number of measurement and the dimension of unknown θ are considered as m, n 301, 301 and the elapsed time T is 15 s.The magnitude of the solution in Figure 3 is unrealistically high and thus this solution is completely useless.Figure 4 shows the distribution of the singular values from the singular value decomposition 28 of the approximated system.Some of singular values become very small, and this makes the original mathematical model very unstable.In order to address this difficulty, we follow the solution procedure by formulating the stochastic inverse problem.
The hierarchical Bayesian formulation 5.3 is adopted as a probabilistic model for the stochastic inverse solution.In order to explore the posterior state, Markov chain Monte Carlo is used as described in Section 6.For the algorithm, the initial guesses λ 0 and σ 0 for hyperparameters are taken to be 1 and 0.1 and θ 0 is considered to be zero vectors.In addition, the pairs of parameters α 1 , β 1 and α 2 , β 2 are taken to be 1, 0.1 and 1, 0.1 , respectively.The total number of samples N MCMC for Algorithm 1 is taken to be 50,000, and the last 25,000 realizations are used to estimate statistics such as mean and standard deviation for the nonlinear damping θ.
The numerical results are shown in Figure 5.The upper and lower dotted lines denote the 95% credible interval.It is easily found that the posterior mean E{θ} is fairly accurate compared with exact solution, and the credible interval shows the quantification of its associated uncertainty corresponding to the measured data.Note here that the number of samples means the number of realization for the unknown θ {θ 1 , θ 2 , . . ., θ m }, and for each element θ i of the unknown θ has 25,000 realizations.Examples of trace plots are illustrated in Figure 6 for θ 1 and θ 100 : θ 1 and θ 100 correspond to the samples for t 0 s and 5 s, respectively.The quality of the MCMC sample may be estimated by inspecting the autocorrelation of a solution sequence 13 .The scale parameter σ q for a proposal distribution can affect the efficiency of the Markov chain.The stronger correlation may result in poor sampling efficiency.Here, we used σ q 0.05 to produce a sample with reasonably low autocorrelation.Now it is time to determine the nonlinear damping.Once the solution is obtained, the nonlinear damping can be determined from the relationship θ t −B y t , ẏ t .In this study, the posterior mean E{θ} is used and the result is depicted in Figure 7. Finally, the displacement and velocity are resimulated by using the identified nonlinear damping in Figure 8.It is confirmed that the results are in well coincidence with exact motion response.

A Particular Application to Realistic Problem: Ship Roll Motion
The results in the numerical example illustrated in Section 7 are fairly good.However, the results have a limitation that the order of damping magnitude is same order with the restoring.In practice, the order of damping magnitude is so low in relative to restoring force; the identification of damping is thus always more problematic.Considering this point, we also applied the present identification procedure to a realistic problem related to ship roll motion by conducting some experiments.The roll motion of ships, barges, and similar floating structures are highly involved in strong nonlinearity.This is mainly because of the complex interaction between floating body and surrounding fluid; see Figure 9.The more important thing is that the dynamic stability of floating body in realistic sea is dependent on its rolling motion.For this purpose, accurate prediction of motion response to various loading conditions is necessary and it requires the precise information on roll damping.

Experimental Set-Up
All experiments were conducted in Ocean Engineering Basin at the University of Tokyo.The basin, which is also called as towing tank, is designed to perform various tests related to various kinds of floating structures.A model tested is shown in Figure 10 and its particulars are summarized in Table 1.For experiments, the model ship was restrained in all degrees of

Analysis of Model Test Data
The mathematical model of ship rolling in calm water is generally written as a second-order differential equation: where I represents the actual mass moment of inertia; A 44 is the added mass moment of inertia; B represents the nonlinear roll damping; C 44 represents the coefficient of restoring moment that is hydrostatic force induced by hydrostatic pressure the static buoyancy force exerted by a fluid at equilibrium due to the force of gravity.The coefficient C 44 is simply expressed as the multiplication of the displacement of ship and the distance between where ω n is the resonance frequency of a ship and it can be obtained by analyzing the free-roll decay curve in Figure 11.For the test model, ω n 6.905 rad/s.In a similar way described in previous section, we first start with the roll response data in Figure 11.The calculated pseudo-displacement g is shown in Figure 12.Now, we can construct the approximate system as in 3.1 .Figure 13 shows the solution through pseudo inverse and it turns out to be very illconditioned.To address this difficulty, we also try to consider the original ill-conditioned system in the stochastic space by following the procedures described in Section 4. By adopting the hierarchical Bayesian formulation 5.3 , the original system can be treated as well-posed system in an expanded stochastic space.The solution θ can be obtained through Markov chain Monte Carlo and illustrated in Figure 14.
The nonlinear damping in Figure 15 can be identified by the relationship θ t − B φ t .Sometimes, it may be convenient to express nonlinear damping as a class of functions rather than a set of data.We caution that it is sometimes extremely difficult to find an analytic function for nonlinear damping as illustrated in the numerical example in Section 7.For such cases, we have to pay attention to choose a class of functions for the identified damping.However, for the present example, we can find that the damping is mainly dependent on the roll angular velocity.In this case, we can approximate the nonlinear damping as a polynomial of the roll angular velocity depicted in Figure 15 through the least squares method which minimizes the residual The solid-line in Figure 15 shows a polynomial approximation of the third order for the identified nonlinear damping: p φ 0.2144 φ3 0.0264 φ2 0.3113 φ.If we know the analytical solution of the nonlinear damping, the order of polynomial can then be decided by comparing the norm between the analytic solution and the polynomial approximation p n φ with respect to the order n.However, the analytical solution to the nonlinear damping B is not available for this application.To ensure the validity of the identified solution, we compare the resimulated roll response using identified nonlinear roll damping with the measured roll response data.For the simulation, the Runge-Kutta method is used.The result is shown in Figure 16.The result is well coincident with the measured roll response.This proves the validity of the present identification procedure and thus the polynomial approximation.

Conclusions
An original output-only and nonparametric procedure for the identification of nonlinear damping in nonlinear oscillatory motion was proposed.The nonlinear damping identification was formulated as a stochastic inverse problem defined on the state space.Probabilistic modeling for the stochastic inverse problem was developed.The way to design computational tools for stochastic inverse solutions with full-probabilistic specification was also illustrated.Numerical experiments were made to show the workability and applicability of the proposed method.In addition, the present procedure was also applied to a realistic problem related to ship roll motion.From the results, it was concluded that the proposed method is well justified for detecting the nonlinear damping in the nonlinear oscillatory system.Through the results, it was also shown that the stochastic inverse formalism has lots of distinct features over deterministic inversion.Although this study has some limitations that the only nonlinear single degree of system is considered, the results of the present work may give an insight applicable to many other nonlinear system identification procedures.

Figure 1 :
Figure 1: The simulated responses of the numerical example.

Figure 2 :Figure 3 :
Figure 2:The calculated pseudo-displacement g and noisy data g δ for the first example.

Figure 4 :
Figure 4: The distribution of the singular values.

Figure 5 :
Figure 5: The numerical results from MCMC computations.

Figure 7 :
Figure 7: a The identified nonlinear damping 3D plots and b its 2D projection on the planes y, B and dy/dt, B .

Figure 12 :Figure 13 :
Figure 12: The calculated pseudo-displacement g for the identification of roll damping of a ship.

Figure 14 :Figure 15 :
Figure 14: The numerical results from MCMC computations for the identification of roll damping of a ship.

Figure 16 :
Figure 16:The resimulated responses using the identified nonlinear damping.

Table 1 :
Particulars of the model ship.