Abstract
The determination of an external force is a very important task for the purpose of control, monitoring, and analysis of damages on structural system. This paper studies a stochastic inverse method that can be used for determining external forces acting on a nonlinear vibrating system. For the purpose of estimation, a stochastic inverse function is formulated to link an unknown external force to an observable quantity. The external force is then estimated from measurements of dynamic responses through the formulated stochastic inverse model. The applicability of the proposed method was verified with numerical examples and laboratory tests concerning the wave-structure interaction problem. The results showed that the proposed method is reliable to estimate the external force acting on a nonlinear system.
1. Introduction
In the field of engineering, it is often very important to estimate external loads acting on dynamic structural systems for a design and analysis of the structural systems. The proper determination of external loads may contribute to extensive practical applications in control, monitoring, and analyzing engineering systems such as bridges under static or dynamic loading, floating structures on waves, and supply systems such as pipes.
However, for some structural systems, it is sometimes difficult to directly measure an external force for some reasons such as installation difficulty of the measurement devices and large magnitude forces. These difficulties led to a necessity of alternative methods to estimate external forces. Therefore, many methods have been developed by means of the inverse problem formalism, inferring an unobservable value from an observable value, except very few cases where the force can be measured directly.
Stevens [1] comprehensively reviewed force identification problems that arise both in laboratory experiments and field applications. Moreover, it was explained that there exist difficulties in obtaining a satisfactory result since the problem is generally unstable. General numerical methods, which proceed in a step-by-step fashion, fail to treat such a problem in a reliable and practical way. To tackle this difficulty, Zhu and Law [2] investigated a method based on the modal superposition and the regularization technique. Doyle [3] and Inoue et al. [4] studied a reconstruction technique based on the least squares method using the singular value decomposition (SVD). Gunawan et al. [5] developed a new method based on the quadratic spline approximation to solve the problem.
The aforementioned studies assume that structural systems of interests are linear. However, in practice, nonlinearities are always present in various forms for mechanical or structural systems. Therefore, it is necessary to consider these nonlinearities for reliable system design and analysis. To cite a few examples of the force identification for nonlinear systems, Ma and Ho [6] developed an inverse method that comprises the extended Kalman filter and the recursive least-squares estimator. Jang et al. [7, 8] presented an inverse procedure based on a deterministic classical regularization method, the Landweber iteration.
It is worth here noting that the aim of the above methods, regardless of considering a linear or nonlinear system, is to make the solution stable by modifying the original formulation. These methods yield only a single estimate of unknowns without quantifying associated uncertainties in the solution. Thus, these kinds of methods are often referred to as the deterministic inverse method.
Due to increasing demands on robustness and reliability as well as due to the rapid growth of computational capacity, it has been more and more important to solve problems under conditions of uncertainty. As a result, recently, problems considering uncertainties have thus attracted much research attention in diverse fields such as electric impedance tomography [9] and heat conduction [10, 11]. There are a few cases concerning the force identification with rigorous consideration of uncertainties. Wu and Law [12] recently developed a moving force identification technique based on a statistical system model. They proved the workability of the technique by applying it to the bridge-vehicle system with the Gaussian system uncertainties arising from road surface roughness. Wu and Law [13] also proposed a dynamic analysis method on the bridge-vehicle interaction problem with the non-Gaussian system uncertainties.
In this paper, an original method based on a stochastic inversion is developed to determine external forces acting on a nonlinear system from measurements of a system response. The method studied here comprises the following four parts: (i) construction of a mathematical inverse model linking an external force to a motion response, (ii) formulation of a stochastic inverse function that allows to determine an external force, (iii) design of probabilistic model to quantify various uncertainties arising from an insufficiency of information on parameters of interests and measurement errors, and (iv) computation of inverse solutions by designing computational tools with full probabilistic description.
In the first part of this paper, the mathematical formulation on the proposed method is presented. Complete description of the proposed method is illustrated through numerical examples. In addition, mathematical characteristics of the constructed model are also discussed. In the second part of this paper, the proposed method is validated through a practical application to a real problem regarding the wave-structure interaction. Experiments for a ship subjected to wave excitation are performed in the ocean engineering basin. Results of the analysis of the experimental data show that wave excitations can be successfully estimated from measurements of the response alone, using the proposed method.
2. Mathematical Formulation
2.1. Preliminaries
The following nonlinear physical model is first considered: where is the mass, is the damping, is the restoring force, and is the nonlinear system including nonlinear damping and restoring forces. Initial conditions for this system are given as Let be basic solutions of the following homogeneous equation associated with (2.1): where and . Then, the solution of (2.1) is given in the form where ’s are continuous functions that satisfy The first and second derivatives of (2.4) are Substituting (2.6) into (2.1) leads to Equation (2.3) implies that the first part of (2.7) is zero since ’s are the solutions of (2.3). Then, (2.7) can be rewritten as Equations (2.5) and (2.8) lead to the following system of equations: Therefore, ’s are given by The inverse matrix in the above equation always exists because ’s are linearly independent. Integrating (2.10) gives where is Wronskian defined as and ’s are arbitrary additive constants. Substituting (2.11) into (2.4) leads to the following equation [14, 15]:
If , are chosen by satisfying , and , , then the constants , become , . The preceding derivation is based on the concept of variation of parameters. The reader can refer to [14, 15] for more detailed information.
In (2.12), can be considered as an effective force acting on the nonlinear vibrating system. Part of the importance of this integral formulation lies in its use not only for estimating the external forces but also for determining the nonlinear system that; characterizes dynamic oscillations. For the present study, the aim is to find the external force acting on nonlinear system. Thus the nonlinear system is assumed to be known.
2.2. A Stochastic Inverse Model for an Indirect Measurement
If a motion response is measured for the time duration and the nonlinear system is known, (2.12) can be rewritten [7, 8] as where and the operator is defined by with the kernel The estimation of the external force can be achieved by solving (2.13) using the measured response data. This problem can be interpreted as an inverse problem of inferring the unknown value from the observable physical quantity .
Equation (2.13) is classified as the integral equation of the first kind [16], whose unknown is only in the integrand. According to the inverse problem theory [16, 17], the first-kind integral equation is typically ill posed in the sense of stability. The approximation of the integral equation of the first kind yields an ill conditioned system regardless of the choice of approximate methods such as the quadrature method and the Galerkin method with orthonormal basis functions. This ill-conditioning causes a loss of accuracy for direct factorization solution methods. Consequently, conventional or standard methods in numerical algebra cannot be used in a straightforward manner for solving numerical system of (2.13). Thus, it is necessary to adopt more sophisticated methods for stable solution procedure.
In this study, we formulate a stochastic inverse function to achieve a stable solution procedure. By considering all variables of interests to be random variables, (2.13) can be recast in the form of a stochastic function. Note that, from now on, the capital letter and lower letter refer to a random variable and its realization, respectively.
For an arbitrary nonempty sample space , let us define the multivariate random variable , , on the state space Then the random variable , , can be defined to be for all possible subsets in on the state space Therefore, for each realization , of , the following relationship is given: In the preceding definitions, is the number of unknowns, is the number of measurements, and is real Euclidean space.
Equation (2.20) can be interpreted as the stochastic function corresponding to (2.13) because it is a function of random variables. Thus, the problem of finding from a single realization is the stochastic inverse problem.
2.3. Probabilistic Modeling of Stochastic Inverse Solution
A single realization , which is the directly observable value, can be readily calculated through (2.14) when the response is measured. Using this value , the unknown of the stochastic inverse problem (2.20) can be expressed [17, 18]: where is the likelihood function, is the prior probability density function, and is the normalizing constant.
The probabilistic expression (2.21), which is known as the posterior probability distribution function (PPDF), can be considered as the solution to the stochastic inverse problem (2.20). However, (2.21) does not give any information on the desired value in its present form. It is necessary to design a probability model according to the proper physical process to have a meaningful probabilistic expression. This modeling enables to quantify uncertainties in the solution associated with insufficient information on the system and measurement errors. From now on, the symbol is used in place of .
For the purpose of modeling, components in the right-hand side of (2.21) are separately considered. At first, it is worth noting that the normalizing constant is not necessary to be computed for sampling procedures [9–11, 17, 18] such as Monte Carlo simulation, which is used to estimate the target probability distribution. The PPDF (2.21) can thus be evaluated with the slightly simplified form:
If a measurement error in the observed data is assumed to be independent, identically distributed Gauss random noise with zero mean and variance , the likelihood function can be modeled in the form where refers to Euclidean norm.
Finally, it needs to consider the prior distribution , reflecting the knowledge on the unknowns before the data is measured. In this study, a pairwise Markov random field [17] is adopted for prior distribution modeling. The MRF model for the unknown is of the form where the matrix is determined by In (2.25), is the number of neighbors for the point , means that the points , are adjacent.
Using the models for the likelihood (2.23) and the prior (2.24), the PPDF (2.21) is given as Here, it is worth discussing the role of parameters , . The scaling parameter controls the prior model. The variance explains the noise level in the measurement data . These parameters play an important role in the solution procedure because the PPDF is characterized by these values as in (2.26). If the variance of the measurement data and the scaling parameter are given, the stochastic inverse solution is readily estimated with the probability model (2.26). However, in practice, it is difficult to know these parameters a priori unless experiments are repeatedly conducted enough to collect the data.
The choice of hyperparameters and sometimes causes another difficulty in solution procedure. However, this difficulty can be naturally resolved by introducing a hierarchical model for the PPDF [10, 11, 19]: with the prior density and .
2.4. Computation of Stochastic Inverse Solution
Equations (2.26) and (2.27) obviously contain the information on the desired solution, that is, the external force . To extract useful information, these probabilistic expressions need to be evaluated through a numerical sampling technique such as Markov chain Monte Carlo (MCMC) simulation [20, 21].
The MCMC methods are a class of algorithms for sampling from probability distributions based on constructing a Markov chain, which is characterized by the rule that the next state depends only on the current state. The generated samples from the algorithm are used to approximate the target distribution. Metropolis-Hastings and Gibbs algorithms [22, 23] are most widely used for sampling from multidimensional distributions, especially when the number of dimensions for unknowns is high. A sequence of random samples from the PPDF (2.27) can thus be obtained by the following hybrid algorithms [22, 23]. (i) Initialize , and (ii) For . Sample Sample Sample Sample Sample if else Sample Sample if else.In the above algorithm, is the total number of samples, is the uniform distribution, and is an easy-to-sample proposal distribution with or .
3. Numerical Investigations
In this section, we want to numerically investigate the proposed method to determine external forces acting on nonlinear system through a particular numerical example.
3.1. Numerical Example
As a numerical example, the following physical model is considered: with the physical values; , and . Comparing (3.1) with (2.1) leads to . For the purpose of testing the proposed method, an external force, which is a repeated nonharmonic forcing, is given by where the values are , and .
At first, the motion responses were simulated with zero initial conditions and by solving (3.1) through the numerical integration scheme such as the fourth-order Runge-Kutta method. Figure 1 shows the exact external force in (3.2) and the corresponding responses. The dynamic responses in Figure 1 were assumed as the measured data for the inverse procedure. For the numerical example, the number of measurements and unknown was taken to be .
(a)
(b)
(c)
The aim is now to estimate the external force from the measurement of system responses using the proposed stochastic estimation method in Section 2. With the measurements in Figure 1, observable quantities were computed through (2.14). To examine the effect of the level of noises in data, two noisy data with the noise level of and were generated by adding independent random noises to the generated responses. Figure 2 shows the exact and two noisy data for .
| (a) |
| (b) |
3.2. Singular Value Analysis
In this subsection, the characteristic of the deterministic inverse model (2.13) is examined by singular value analysis [24], which provides quantitative measures for assessing the numerical system. For this purpose, the discretized or approximated form of (2.13) can be written as or simply The matrix is thus expressed as a lower triangular Toeplitz matrix: By using the singular system [16, 17, 24], the matrix can be decomposed as where is an unitary matrix, the matrix is an diagonal matrix with nonnegative real numbers , which are known as the singular values of , and is the conjugate transpose of . With this singular system, the inverse solution can easily be written in the form
Equation (3.7) is known as the pseudo inverse [24], which is also referred to as a generalized inverse. The purpose of constructing a pseudo inverse in (3.7) is to obtain a systematic expression that can serve as the inverse in some sense for a wider class of numerical systems than invertible ones. The inverse solution was computed through the pseudo inverse (3.7) with the observable quantity . The results are shown in Figure 3. It can be seen that the computed solutions are totally unstable and thus useless.
| (a) |
| (b) |
The cause and the degree of instability in the solution can be understood by analyzing the behavior of singular systems, more specifically, each component in the right-hand side of (3.7). Figure 4 shows behaviors of , , and the ratio . It can be found that the coefficients become larger than most of singular values . Thus, the inverse solution , which is the sum of the ratio , is dominated by the small singular values. This explains the reason for significant sign changes in Figure 3. The instability, frequent sign changes, may be accelerated and amplified with the increasing amount of noise in data. This can be easily checked by comparing Figures 3(a) and 3(b). The magnitude of the solution in Figure 3(b) is larger than the one in Figure 3(a).
| (a) |
| (b) |
In summary, for a stable solution, the approximated system should satisfy the condition that coefficients should decay to zero faster than singular values . This condition is known as the discrete Picard condition [24]. The figure describing the behavior of singular values is also known as the Picard plots. Since the deterministic inverse model (2.13) fails to satisfy this condition, it is ill posed in the sense of stability. Therefore, the conventional or standard methods in numerical algebra cannot be used to solve the inverse model (2.13).
3.3. Solution to Stochastic Inverse Model
As pointed out, the inverse model (2.13) loses its stability when it is approximated. Thus, the usual numerical scheme cannot be used in straightforward manner. For stable solution procedure, the stochastic inverse method described in Section 2 is followed here by transforming the original deterministic inverse model into stochastic state space.
Accordingly, the inverse solution, unknown external force , has the form as the PPDF in (2.27). To extract necessary information from the PPDF, the MCMC algorithm described in Section 2.4 is used. For the algorithm, the values , and were initialized as , , and zero vectors for , respectively. The pairs of parameters and in (2.27) determine the initial shape of prior densities and . In the MCMC algorithm, these prior densities are greatly refined through the conditioning on the data regardless of the initial choice. In this study, the pairs of parameters and were taken to be and to impose nonnegativity. The total number of samples for the algorithm was taken to be , and the last realizations were used to estimate statistics such as mean and standard deviation for external force since it allows the beginning of the Markov chain. This beginning of chain is often called the burn-in. The scale parameter for a proposal distribution can affect the efficiency of the Markov chain. In this study, was used to produce a sample with reasonably low autocorrelation.
Numerical results are shown in Figure 5. The upper and lower dotted lines denote the 95% credible interval. The credible interval quantifies the degree of uncertainties in the solution. Comparing Figures 5(a) and 5(b) implies that the degree of uncertainties depends on the noise level in the measurement data. The credible interval in Figure 5(a) is narrower than the one in Figure 5(b). It can also be seen that the posterior mean is fairly accurate compared with the exact solution, regardless of the level of noise.
| (a) |
| (b) |
Figure 6 shows examples of trace plot with respect to the number of sampling for the case . These points and are highlighted with arrows in Figure 5(a). One of the features that can be observed in Markov chain is that it takes a while to properly sample the target distribution. The evolution of components may be a good tool to check if the burn-in begins. The curve in Figure 6 clearly shows the burn-in phase. This implies that the chain has reached its stationary state and thus works properly.
(a)
(b)
(c)
(d)
4. A Particular Application: Estimation of Wave Exciting Forces Acting on a Ship
The results in the previous section were fairly good. However, it has a limitation that numerical simulations were conducted based on the synthesized data. In this section, to ensure the practicability, the proposed method is also applied to the real practical problem regarding the problem of estimating wave exciting forces on a ship in roll motion. The estimation of wave exciting forces is one of the most important issues in offshore structure design since waves are the most important environmental phenomena that affect the system stability of offshore structures.
4.1. Motion Equation of a Ship in Roll Motion
The mathematical model of ship rolling in waves (Figure 7) is generally expressed as the following nonlinear differential equation: where is the actual mass moment of inertia, is the added mass moment of inertia, is the nonlinear roll damping function, is the coefficient of restoring moment, is the roll response, and represents a wave exciting moment acting on a ship. The sum of and is often referred to as the virtual mass.
Equation (4.1) looks similar to the general nonlinear physical model (2.1). However, the mechanism behind this is somewhat complicated because the interaction between the ship hull and the surrounding water should be taken into account. This interaction causes the mechanism for damping and restoring forces in ship motions. Figure 8 illustrates an analogy of the mechanical system and the motion of a ship in water.
(a)
(b)
More specifically, forces acting on ships can be decomposed into hydrostatic and hydrodynamic forces. Hydrostatic force is the force induced by hydrostatic pressure (the static buoyancy force) exerted by a fluid at equilibrium due to the force of gravity. This hydrostatic force is related to the restoring mechanism of a ship in water.
Hydrodynamic forces can be considered as the forces resulted from the pressure field disturbed by a moving body. Hydrodynamic forces are divided into three main components. The first is the Froude-Krylov force due to the pressure field in the incident wave. The second is the diffraction force due to the scattering of the incident wave field. The last is the radiation force due to the waves generated by the oscillatory motions of the body. Among them the radiation force is related to the added mass and damping mechanism of a ship motion. The wave exciting force is concerned with the Froude-Krylov force and the diffraction force. For more detailed information, the reader can refer to [25].
In this section, the attention is focused on the estimation of the wave excitation . For an application, all physical values are separately identified a priori except the wave excitation.
4.2. Experimental Setup
For the purpose of practical application, laboratory tests regarding ship rolling motions in waves were conducted in the Ocean Engineering Basin at The University of Tokyo. The basin, which is often called the towing tank, is designed to perform tests related to various kinds of floating structures. Figure 9 shows the test model and its body plan illustrating the shape of hull form. Table 1 summarizes particulars of the test model.
(a)
(b)
Figure 10 shows an overview of experimental setups. Regular wave fields with a fixed frequency were generated by using the plunge type wave-maker. Then motion responses were measured with the potentiometer attached to the center of gravity of the test model. The incident wave was also recorded by far and near wave probes from the test model. During experiments, the test model was located in the middle of the basin to minimize the effects of reflected waves. For simplicity, other motions of the test model were restricted except the roll motion.
Figure 11 shows a typical example of the measured data. We used the stationary region of the data highlighted with red box in Figure 11 to estimate wave exciting moment. For demonstration purposes, three frequency regions were chosen-low, high, and near resonance frequency.
(a)
(b)
(c)
4.3. Analysis of the Experimental Data
As a first attempt, the proposed method was applied to the roll response data when the frequency in Figure 12. The data was measured with the sampling period s. Using this data, the observable quantity in Figure 13 was calculated through (2.14). This quantity leads to the construction of the PPDF for wave exciting moment as in (2.27).
In the same way used in the previous section, the information on unknown wave exciting moment was extracted via MCMC algorithm. Without loss of generality, the values , and were initialized as , , and zero vectors for . The pairs of parameters and were also taken to be and .
Figures 14 and 15 show, respectively, the solution obtained by the pseudo inverse (3.7) and the proposed stochastic inverse method. It can be seen that the stochastic inverse solution is stable while the pseudo inverse solution loses its stability and accuracy. Moreover, the result in Figure 14 shows a clear periodicity as expected from the periodic excitation. It is also observed that the stochastic inverse solution provides a credible interval explaining uncertainties in unknown quantity .
For the purpose of validation, roll response was resimulated using the posterior mean of the identified wave exciting moment since the exact solution for cannot be known a priori for the real practical problems unlike numerical simulations. Figure 16 shows the comparison between the re-simulated and the measured responses. The result shows that the re-simulated response is well coincident with the measured one. This proves the validity of the estimated wave exciting moment for this case.
The same procedure was also applied to the other response data when and . All results are shown in Figures 17 and 18, respectively. In both cases, it can be seen that the re-simulated responses using the estimated wave exciting moment well agree with those of measurements.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
When calculating the wave exciting force acting on a ship, the strip theory [25] is widely used in the field of naval architecture and ocean engineering. The strip theory is a popular approximation tool to calculate hydrodynamic forces for ships or marine vehicle. The key idea behind this theory is that hydrodynamic forces can be estimated based on a stripwise computation by slicing the entire body of a ship in the longitudinal direction. Resultant wave force is then expressed as the sum of loads that is calculated in each segment on the sliced body.
To ensure the validity of the present method, a qualitative comparison is also made in Figure 19 for wave exciting moments estimated by three different sources: the proposed method, the strip theory, and the experiments. The horizontal axis represents the ratio of the length of test model and wave length. The vertical axis represents nondimensional wave exciting moment. As can be seen in Figure 19, all the results well agree with each other. However, it can be found that there exist slight differences between the proposed method and that of the strip theory when is small or large. This may be mainly because of the linear assumption in the strip theory.
4.4. Irregular Wave Excitation
It is worth emphasizing that the proposed method is not limited to the case of the stable sinusoidal responses. It is also applicable to the estimation for any form of wave excitation, provided that the motion response is properly recorded. To test this, we also conducted another experiment regarding the ship motion subjected to the wave excitation induced by irregular waves.
Figure 20 shows a time history of measurements. It can be observed that the response is not simple harmonic since the incident wave is irregular. This roll response was then used, according to the proposed stochastic inverse procedure, to estimate the irregular wave excitation. Figure 21 shows the estimated results. Finally, using this estimated wave excitation, the motion was re-simulated numerically in the same manner described in the previous examples. The result in Figure 22 well coincides with the measurement. This proves the applicability of the proposed method to other types of waveforms.
(a)
(b)
(c)
5. Conclusions
A new method has been proposed to estimate the external forces acting on a nonlinear vibrating system, based on formulating the stochastic inverse model from the measured roll response data. The method has been evaluated through an analysis of some digitally simulated data for a numerical example. It has been shown that the proposed method can be used to not only identify the external force but also detect uncertainties in solution arising from measurement error or insufficient information on the system.
To ensure applicability, the wave-structure interaction problem regarding the roll motion of a ship has also been treated as a practical application. To this end, a series of experiments were conducted in the ocean engineering basin. The analysis results demonstrate that the proposed method has been successfully applied to estimate the various types of wave excitation. The estimated results were qualitatively and quantitatively good in all cases.
Acknowledgments
The authors would like to thank the Editor and anonymous reviewers for their valuable comments, suggestions, and constructive criticism, which were very helpful in improving the paper substantially.