Abstract

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.

1. 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 in-service 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.

2. Mathematical Description

Consider nonlinear single degree of freedom system as
with initial conditions
In (2.1), and are nonlinear damping and restoring force, respectively. In this study, the nonlinear restoring function is considered to be known a priori and assumed to have the form
where and are linear and nonlinear restoring components. Then, (2.1) can be rewritten as
From the concept of variation of parameters, nonlinear motion equation (2.1) can be transformed to following nonlinear Volterra integral equation of the second kind [21]:
where and satisfy
and Wronskian is

3. 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 . Let denote the measured response data, that is, ; then we can obtain the following system for unknown nonlinear damping :
where

The identification of the unknown nonlinear damping 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.

4. 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 and in (3.1) as random variables. Let us define a multivariate random variables ; then a multivariate random variable can be defined by
For each realization of , the following relationship is given:
where 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 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 is observed, then a single realization of directly observable parameter can be obtained from (3.3), and then the solution to the stochastic inverse problem can be expressed by Bayes’ rule [11–14]:
where , and 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 in (4.3) can then be evaluated as follows:

It is usually easy to obtain the likelihood 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 and refer to Euclidean norm and the number of measurements, respectively.

Next, we should consider the prior probability density function 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, is dimension of , is the scaling parameter, and the matrix is determined as
where is the number of neighbors for the point . Then, with the likelihood (4.5) and the prior distribution (4.6), posterior probability density function (4.3) can be formulated as

5. 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 . 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
where and are pairs for gamma distribution for and inverse gamma distribution for :

6. 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, is the total number of samples, refers to , is an easy-to-sample proposal distribution, and the conditional which follows a multivariate Gaussian distribution can be expressed by
In a similar manner, the full conditional and can also be derived as

Algorithm 1

7. 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 . Comparison of (2.1) and (7.1) shows that and .

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 and . 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 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 with 10.1% noise is also plotted in Figure 2.

Figure 1: The simulated responses of the numerical example.

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

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 and the elapsed time 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.

Figure 3: The numerical solution of the numerical example through pseudoinversion.

Figure 4: The distribution of the singular values.

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 and for hyperparameters are taken to be 1 and 0.1 and is considered to be zero vectors. In addition, the pairs of parameters and are taken to be and , respectively. The total number of samples 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 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 , and for each element of the unknown has 25,000 realizations. Examples of trace plots are illustrated in Figure 6 for and : and correspond to the samples for and , respectively. The quality of the MCMC sample may be estimated by inspecting the autocorrelation of a solution sequence [13]. The scale parameter for a proposal distribution can affect the efficiency of the Markov chain. The stronger correlation may result in poor sampling efficiency. Here, we used to produce a sample with reasonably low autocorrelation.

Figure 5: The numerical results from MCMC computations.

Figure 6: The trace plots for and and their posterior distributions.

Now it is time to determine the nonlinear damping. Once the solution is obtained, the nonlinear damping can be determined from the relationship . In this study, the posterior mean 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.

Figure 7: (a) The identified nonlinear damping 3D plots and (b) its 2D projection on the planes and .

Figure 8: The resimulated responses using the identified nonlinear damping in Figure 7.

8. 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.

Figure 9: Interaction between ship hull and surrounding water.

8.1. 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 freedom during experiments, excepting roll motion. The model ship is placed at the center of the basin and roll angle is measured with potentiometer attached to the center of gravity of the model ship.

Table 1: Particulars of the model ship.

Figure 10: Test model of vessel.

8.2. Analysis of Model Test Data

The mathematical model of ship rolling in calm water is generally written as a second-order differential equation:
where represents the actual mass moment of inertia; is the added mass moment of inertia; represents the nonlinear roll damping; 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 is simply expressed as the multiplication of the displacement of ship and the distance between metacentric height and center of gravity GM: . By eliminating the inertia coefficient, (8.1) can be rewritten as
where 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, rad/s.

Figure 11: Roll responses of free-decay test.

In a similar way described in previous section, we first start with the roll response data in Figure 11. The calculated pseudo-displacement 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 .

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

Figure 13: The numerical solution through pseudoinversion for the identification of roll damping of a ship.

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

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: . 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 with respect to the order . However, the analytical solution to the nonlinear damping 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.

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

9. 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.