# Investigation of a Stochastic Inverse Method to Estimate an External Force: Applications to a Wave-Structure Interaction

Authors: S. L. Han; Takeshi Kinoshita
Journal: Mathematical Problems in Engineering
 (2012)

---

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


---

## Body

## 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:(2.1)mx¨(t)+cx˙(t)+kx(t)+SNL(x˙,x;t)=f(t),
where m is the mass, c is the damping, k is the restoring force, and SNL is the nonlinear system including nonlinear damping and restoring forces. Initial conditions for this system are given as
(2.2)x(0)=α,x˙(0)=β.
Let x1,x2 be basic solutions of the following homogeneous equation associated with (2.1):
(2.3)x¨+2ξωnx˙+ωn2x=0,
where ξ=c/2ωn and ωn=k/m. Then, the solution of (2.1) is given in the form
(2.4)x=∑i=12wi(t)xi(t),
where wi’s are continuous functions that satisfy
(2.5)∑i=12w˙ixi=0.
The first and second derivatives of (2.4) are
(2.6)x˙=∑i=12wix˙i,x¨=∑i=12w˙ix˙i+wix¨i.
Substituting (2.6) into (2.1) leads to
(2.7)∑i=12wi(x¨i+2ξωnx˙i+ωn2xi)+∑i=12w˙ix˙i=f-SNL(x,x˙)m.
Equation (2.3) implies that the first part of (2.7) is zero since xi’s are the solutions of (2.3). Then, (2.7) can be rewritten as
(2.8)∑i=12w˙ix˙i=f-SNL(x,x˙)m.
Equations (2.5) and (2.8) lead to the following system of equations:
(2.9)[x1x2x˙1x˙2](w˙1w˙2)=(0f-SNLm).
Therefore, w˙i’s are given by
(2.10)(w˙1w˙2)=[x1x2x˙1x˙2]-1(0f-SNLm)=1W[x˙2-x2-x˙1x1](0f-SNLm).
The inverse matrix in the above equation always exists because xi’s are linearly independent. Integrating (2.10) gives
(2.11)w1(t)=-∫0tx2(τ)mW[f(τ)-SNL(x,x˙)]dτ+D1,w2(t)=-∫0tx1(τ)mW[f(τ)-SNL(x,x˙)]dτ+D2,
where W is Wronskian defined as W(x1,x2)=|x1x˙2-x2x˙1| and Di’s are arbitrary additive constants. Substituting (2.11) into (2.4) leads to the following equation [14, 15]:
(2.12)x=D1x1+D2x2+∫0tx1(τ)x2(t)-x1(t)x2(τ)mW[f(τ)-SNL(x,x˙)]dτ.Ifx1, x2 are chosen by satisfying x1(0)=1, x1˙(0)=0 and x2(0)=0, x2˙(0)=1, then the constants D1, D2 become D1=α, D2=β. 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), f-SNL 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 f but also for determining the nonlinear system SNL that; characterizes dynamic oscillations. For the present study, the aim is to find the external force f acting on nonlinear system. Thus the nonlinear system SNL is assumed to be known.


### 2.2. A Stochastic Inverse Model for an Indirect Measurement

If a motion responsex is measured for the time duration 0<t<T and the nonlinear system SNL is known, (2.12) can be rewritten [7, 8] as
(2.13)g=F(f),
where
(2.14)g(t)=x-αx1-βx2+∫0TK(t,τ)SNL(x,x˙)dτ,
and the operator ℱ is defined by
(2.15)F(*)≡∫0TK(t,τ){*}dτ,
with the kernel
(2.16)K(t,τ)=x1(τ)x2(t)-x1(t)x2(τ)mW.
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 f from the observable physical quantity g.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 F:Ω→RN, F={F1,F2,…,FN}, on the state space
(2.17)ΓF=∏i=1NΓFi⊂RN.
Then the random variable G:Ω→RM, G={G1,G2,…,GM}, can be defined to be
(2.18)G=F(F)
for all possible subsets in Ω on the state space
(2.19)ΓG=∏i=1MΓGi⊂RM.
Therefore, for each realization g={g1,g2,…,gM}, of g∈ΓG, the following relationship is given:
(2.20)g=F(f).
In the preceding definitions, N is the number of unknowns, M is the number of measurements, and R 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 f from a single realization g is the stochastic inverse problem.


### 2.3. Probabilistic Modeling of Stochastic Inverse Solution

A single realizationg, which is the directly observable value, can be readily calculated through (2.14) when the response x is measured. Using this value gobserved, the unknown f of the stochastic inverse problem (2.20) can be expressed [17, 18]:
(2.21)p(f∣gobserved)=p(gobserved∣f)p(f)p(gobserved),
where p(gobserved∣f) is the likelihood function, p(f) is the prior probability density function, and p(gobserved) 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 f 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 g is used in place of gobserved:g=gobserved.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 p(g) 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:
(2.22)p(f∣g)∝p(g∣f)p(f).If a measurement errore in the observed data g is assumed to be independent, identically distributed Gauss random noise with zero mean and variance ve, the likelihood function p(g∣f) can be modeled in the form
(2.23)p(g∣f)∝(ve)-M/2exp⁡(-‖F(f)-g‖222ve),
where ∥·∥2 refers to Euclidean norm.Finally, it needs to consider the prior distributionp(f), reflecting the knowledge on the unknowns before the data g is measured. In this study, a pairwise Markov random field [17] is adopted for prior distribution modeling. The MRF model for the unknown f is of the form
(2.24)p(f)∝γN/2exp⁡⁡(-12γfTWf),
where the matrix W∈RN×N is determined by
(2.25)Wij={nii=j,-1i~j,0otherwise.
In (2.25), ni is the number of neighbors for the point i, i~j means that the points i, j are adjacent.Using the models for the likelihood (2.23) and the prior (2.24), the PPDF (2.21) is given as
(2.26)p(f∣g)∝(ve)-M/2exp⁡(-‖F(f)-g‖222ve)γN/2exp⁡⁡(-12γfTWf).
Here, it is worth discussing the role of parameters ve, γ. The scaling parameter γ controls the prior model. The variance ve explains the noise level in the measurement data g. 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 ve 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 hyperparametersve 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]:
(2.27)p(f,γ,ve∣g)∝(ve)-M/2exp⁡(-‖F(f)-g‖222ve)γN/2exp⁡⁡(-12γfTWf)⋅p(γ)⋅p(ve)
with the prior density p(γ)∝γα1-1e-β1γ and p(ve)∝(ve)-α2-1e-β2/ve.


### 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 f. 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)
Initializef0={f10,f20,…,fN0}, γ0 and ve0(ii)
Fori=0:NMCMC-1.Samplef1i+1~p(f1∣f2i,f3i,…,fNi,γi,vei)Samplef2i+1~p(f2∣f1i+1,f3i,…,fNi,γi,vei)⋮SamplefNi+1~p(fN∣f1i+1,f2i+1,…,fN-1i+1,γi,vei)Sampleu1~U(0,1)Sampleγ*~q(γ*∣γi)ifu1<min{1,p(γ*∣fi+1,vei)q(γi∣γ*)p(γi∣fi+1,vei)q(γ*∣γi)}
γi+1=γ*else
γi+1=γiSampleu2~U(0,1)Sampleve*~q(ve*∣vei)ifu2<min{1,p(ve*∣fi+1,γi)q(vei∣ve*)p(vei∣fi+1,γi)q(ve*∣vei)}
vei+1=ve*else
vei+1=vei.In the above algorithm, NMCMC is the total number of samples, U is the uniform distribution, and q(X*∣Xi) is an easy-to-sample proposal distribution with X=γ or ve.


## 2.1. Preliminaries

The following nonlinear physical model is first considered:(2.1)mx¨(t)+cx˙(t)+kx(t)+SNL(x˙,x;t)=f(t),
where m is the mass, c is the damping, k is the restoring force, and SNL is the nonlinear system including nonlinear damping and restoring forces. Initial conditions for this system are given as
(2.2)x(0)=α,x˙(0)=β.
Let x1,x2 be basic solutions of the following homogeneous equation associated with (2.1):
(2.3)x¨+2ξωnx˙+ωn2x=0,
where ξ=c/2ωn and ωn=k/m. Then, the solution of (2.1) is given in the form
(2.4)x=∑i=12wi(t)xi(t),
where wi’s are continuous functions that satisfy
(2.5)∑i=12w˙ixi=0.
The first and second derivatives of (2.4) are
(2.6)x˙=∑i=12wix˙i,x¨=∑i=12w˙ix˙i+wix¨i.
Substituting (2.6) into (2.1) leads to
(2.7)∑i=12wi(x¨i+2ξωnx˙i+ωn2xi)+∑i=12w˙ix˙i=f-SNL(x,x˙)m.
Equation (2.3) implies that the first part of (2.7) is zero since xi’s are the solutions of (2.3). Then, (2.7) can be rewritten as
(2.8)∑i=12w˙ix˙i=f-SNL(x,x˙)m.
Equations (2.5) and (2.8) lead to the following system of equations:
(2.9)[x1x2x˙1x˙2](w˙1w˙2)=(0f-SNLm).
Therefore, w˙i’s are given by
(2.10)(w˙1w˙2)=[x1x2x˙1x˙2]-1(0f-SNLm)=1W[x˙2-x2-x˙1x1](0f-SNLm).
The inverse matrix in the above equation always exists because xi’s are linearly independent. Integrating (2.10) gives
(2.11)w1(t)=-∫0tx2(τ)mW[f(τ)-SNL(x,x˙)]dτ+D1,w2(t)=-∫0tx1(τ)mW[f(τ)-SNL(x,x˙)]dτ+D2,
where W is Wronskian defined as W(x1,x2)=|x1x˙2-x2x˙1| and Di’s are arbitrary additive constants. Substituting (2.11) into (2.4) leads to the following equation [14, 15]:
(2.12)x=D1x1+D2x2+∫0tx1(τ)x2(t)-x1(t)x2(τ)mW[f(τ)-SNL(x,x˙)]dτ.Ifx1, x2 are chosen by satisfying x1(0)=1, x1˙(0)=0 and x2(0)=0, x2˙(0)=1, then the constants D1, D2 become D1=α, D2=β. 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), f-SNL 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 f but also for determining the nonlinear system SNL that; characterizes dynamic oscillations. For the present study, the aim is to find the external force f acting on nonlinear system. Thus the nonlinear system SNL is assumed to be known.


## 2.2. A Stochastic Inverse Model for an Indirect Measurement

If a motion responsex is measured for the time duration 0<t<T and the nonlinear system SNL is known, (2.12) can be rewritten [7, 8] as
(2.13)g=F(f),
where
(2.14)g(t)=x-αx1-βx2+∫0TK(t,τ)SNL(x,x˙)dτ,
and the operator ℱ is defined by
(2.15)F(*)≡∫0TK(t,τ){*}dτ,
with the kernel
(2.16)K(t,τ)=x1(τ)x2(t)-x1(t)x2(τ)mW.
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 f from the observable physical quantity g.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 F:Ω→RN, F={F1,F2,…,FN}, on the state space
(2.17)ΓF=∏i=1NΓFi⊂RN.
Then the random variable G:Ω→RM, G={G1,G2,…,GM}, can be defined to be
(2.18)G=F(F)
for all possible subsets in Ω on the state space
(2.19)ΓG=∏i=1MΓGi⊂RM.
Therefore, for each realization g={g1,g2,…,gM}, of g∈ΓG, the following relationship is given:
(2.20)g=F(f).
In the preceding definitions, N is the number of unknowns, M is the number of measurements, and R 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 f from a single realization g is the stochastic inverse problem.


## 2.3. Probabilistic Modeling of Stochastic Inverse Solution

A single realizationg, which is the directly observable value, can be readily calculated through (2.14) when the response x is measured. Using this value gobserved, the unknown f of the stochastic inverse problem (2.20) can be expressed [17, 18]:
(2.21)p(f∣gobserved)=p(gobserved∣f)p(f)p(gobserved),
where p(gobserved∣f) is the likelihood function, p(f) is the prior probability density function, and p(gobserved) 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 f 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 g is used in place of gobserved:g=gobserved.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 p(g) 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:
(2.22)p(f∣g)∝p(g∣f)p(f).If a measurement errore in the observed data g is assumed to be independent, identically distributed Gauss random noise with zero mean and variance ve, the likelihood function p(g∣f) can be modeled in the form
(2.23)p(g∣f)∝(ve)-M/2exp⁡(-‖F(f)-g‖222ve),
where ∥·∥2 refers to Euclidean norm.Finally, it needs to consider the prior distributionp(f), reflecting the knowledge on the unknowns before the data g is measured. In this study, a pairwise Markov random field [17] is adopted for prior distribution modeling. The MRF model for the unknown f is of the form
(2.24)p(f)∝γN/2exp⁡⁡(-12γfTWf),
where the matrix W∈RN×N is determined by
(2.25)Wij={nii=j,-1i~j,0otherwise.
In (2.25), ni is the number of neighbors for the point i, i~j means that the points i, j are adjacent.Using the models for the likelihood (2.23) and the prior (2.24), the PPDF (2.21) is given as
(2.26)p(f∣g)∝(ve)-M/2exp⁡(-‖F(f)-g‖222ve)γN/2exp⁡⁡(-12γfTWf).
Here, it is worth discussing the role of parameters ve, γ. The scaling parameter γ controls the prior model. The variance ve explains the noise level in the measurement data g. 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 ve 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 hyperparametersve 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]:
(2.27)p(f,γ,ve∣g)∝(ve)-M/2exp⁡(-‖F(f)-g‖222ve)γN/2exp⁡⁡(-12γfTWf)⋅p(γ)⋅p(ve)
with the prior density p(γ)∝γα1-1e-β1γ and p(ve)∝(ve)-α2-1e-β2/ve.


## 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 f. 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)
Initializef0={f10,f20,…,fN0}, γ0 and ve0(ii)
Fori=0:NMCMC-1.Samplef1i+1~p(f1∣f2i,f3i,…,fNi,γi,vei)Samplef2i+1~p(f2∣f1i+1,f3i,…,fNi,γi,vei)⋮SamplefNi+1~p(fN∣f1i+1,f2i+1,…,fN-1i+1,γi,vei)Sampleu1~U(0,1)Sampleγ*~q(γ*∣γi)ifu1<min{1,p(γ*∣fi+1,vei)q(γi∣γ*)p(γi∣fi+1,vei)q(γ*∣γi)}
γi+1=γ*else
γi+1=γiSampleu2~U(0,1)Sampleve*~q(ve*∣vei)ifu2<min{1,p(ve*∣fi+1,γi)q(vei∣ve*)p(vei∣fi+1,γi)q(ve*∣vei)}
vei+1=ve*else
vei+1=vei.In the above algorithm, NMCMC is the total number of samples, U is the uniform distribution, and q(X*∣Xi) is an easy-to-sample proposal distribution with X=γ or ve.


## 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:(3.1)mx¨+cLx˙+cNLx˙|x˙|+kLx+kNLx3=f(t),
with the physical values; m=1,cL=1,cNL=1,kL=10, and kNL=10. Comparing (3.1) with (2.1) leads to SNL(x˙,x;t)=cNLx˙|x˙|+kNLx3. For the purpose of testing the proposed method, an external force, which is a repeated nonharmonic forcing, is given by
(3.2)f(t)=acos⁡⁡(ωt)+κa23-ϑ24ϑ3cos⁡⁡(2ωt),
where the values are κ=0.1,a=2,ϑ=0.7616, and ω=0.9039.At first, the motion responses were simulated with zero initial conditionsx(0)=0 and x˙(0)=0 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 M and unknown N was taken to be (201,201).External force and corresponding responses.
(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 Section2. With the measurements in Figure 1, observable quantities g were computed through (2.14). To examine the effect of the level of noises in data, two noisy data with the noise level of σ=0.0158 and σ=0.0316 were generated by adding independent random noises e~N(0,σ2) to the generated responses. Figure 2 shows the exact and two noisy data for g.Exact and two noisy data forg.
(a)
σ
=
0.0158
(b)
σ
=
0.0316


### 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
(3.3)g(tj)=∑iΔτK(tj,τi)f(τi),K(tj,τi)=0whentj≤τi,
or simply
(3.4)g=Kf.
The matrix K is thus expressed as a lower triangular Toeplitz matrix:
(3.5)K=[k1k2k1⋮kNkN-1⋱⋯k1].
By using the singular system [16, 17, 24], the matrix K can be decomposed as
(3.6)K=UΣVT=∑i=1NuiσiviT,
where U=(u1,u2,…,uN) is an unitary matrix, the matrix Σ is an diagonal matrix with nonnegative real numbers σi, which are known as the singular values of K, and VT is the conjugate transpose of V=(v1,v2,…,vN). With this singular system, the inverse solution f can easily be written in the form
(3.7)f=∑i=1NuiTgσivi.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 g. The results are shown in Figure 3. It can be seen that the computed solutions are totally unstable and thus useless.Numerical solutions via the pseudo inverse.
(a)
σ
=
0.0158
(b)
σ
=
0.0316The 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 |uiTg|, σi, and the ratio |uiTg|/σi. It can be found that the coefficients |uiTg| become larger than most of singular values σi. Thus, the inverse solution f, which is the sum of the ratio uiTg/σi, 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).Picard plots for the singular system with two noisy data in Figure2.
(a)
σ
=
0.0158
(b)
σ
=
0.0316In summary, for a stable solution, the approximated system should satisfy the condition that coefficients|uiTg| should decay to zero faster than singular values σi. 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 forcef, has the form as the PPDF p(f,γ,ve∣g) 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 γ,ve, and f were initialized as γ0=1, ve0=0.1, and zero vectors for f0, respectively. The pairs of parameters (α1,β1) and (α2,β2) in (2.27) determine the initial shape of prior densities p(γ) and p(ve). 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 (α1,β1) and (α2,β2) were taken to be (1,10) and (1,10) to impose nonnegativity. The total number of samples NMCMC for the algorithm was taken to be 50,000, and the last 25,000 realizations were used to estimate statistics such as mean and standard deviation for external force f since it allows the beginning of the Markov chain. This beginning of chain is often called the burn-in. The scale parameter σX for a proposal distribution q(X*∣Xi)∝exp(-(Xi-X*)2/2σX2) can affect the efficiency of the Markov chain. In this study, σX=0.05 was used to produce a sample with reasonably low autocorrelation.Numerical results are shown in Figure5. 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 E{f} is fairly accurate compared with the exact solution, regardless of the level of noise.Numerical solutions via the proposed method.
(a)
σ
=
0.0158
(b)
σ
=
0.0316Figure6 shows examples of trace plot with respect to the number of sampling for the case σ=0.0158. These points f40 and f140 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.Examples of trace plots for the MCMC simulation and their posterior distributions.
(a)
(b)
(c)
(d)


## 3.1. Numerical Example

As a numerical example, the following physical model is considered:(3.1)mx¨+cLx˙+cNLx˙|x˙|+kLx+kNLx3=f(t),
with the physical values; m=1,cL=1,cNL=1,kL=10, and kNL=10. Comparing (3.1) with (2.1) leads to SNL(x˙,x;t)=cNLx˙|x˙|+kNLx3. For the purpose of testing the proposed method, an external force, which is a repeated nonharmonic forcing, is given by
(3.2)f(t)=acos⁡⁡(ωt)+κa23-ϑ24ϑ3cos⁡⁡(2ωt),
where the values are κ=0.1,a=2,ϑ=0.7616, and ω=0.9039.At first, the motion responses were simulated with zero initial conditionsx(0)=0 and x˙(0)=0 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 M and unknown N was taken to be (201,201).External force and corresponding responses.
(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 Section2. With the measurements in Figure 1, observable quantities g were computed through (2.14). To examine the effect of the level of noises in data, two noisy data with the noise level of σ=0.0158 and σ=0.0316 were generated by adding independent random noises e~N(0,σ2) to the generated responses. Figure 2 shows the exact and two noisy data for g.Exact and two noisy data forg.
(a)
σ
=
0.0158
(b)
σ
=
0.0316


## 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
(3.3)g(tj)=∑iΔτK(tj,τi)f(τi),K(tj,τi)=0whentj≤τi,
or simply
(3.4)g=Kf.
The matrix K is thus expressed as a lower triangular Toeplitz matrix:
(3.5)K=[k1k2k1⋮kNkN-1⋱⋯k1].
By using the singular system [16, 17, 24], the matrix K can be decomposed as
(3.6)K=UΣVT=∑i=1NuiσiviT,
where U=(u1,u2,…,uN) is an unitary matrix, the matrix Σ is an diagonal matrix with nonnegative real numbers σi, which are known as the singular values of K, and VT is the conjugate transpose of V=(v1,v2,…,vN). With this singular system, the inverse solution f can easily be written in the form
(3.7)f=∑i=1NuiTgσivi.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 g. The results are shown in Figure 3. It can be seen that the computed solutions are totally unstable and thus useless.Numerical solutions via the pseudo inverse.
(a)
σ
=
0.0158
(b)
σ
=
0.0316The 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 |uiTg|, σi, and the ratio |uiTg|/σi. It can be found that the coefficients |uiTg| become larger than most of singular values σi. Thus, the inverse solution f, which is the sum of the ratio uiTg/σi, 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).Picard plots for the singular system with two noisy data in Figure2.
(a)
σ
=
0.0158
(b)
σ
=
0.0316In summary, for a stable solution, the approximated system should satisfy the condition that coefficients|uiTg| should decay to zero faster than singular values σi. 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 forcef, has the form as the PPDF p(f,γ,ve∣g) 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 γ,ve, and f were initialized as γ0=1, ve0=0.1, and zero vectors for f0, respectively. The pairs of parameters (α1,β1) and (α2,β2) in (2.27) determine the initial shape of prior densities p(γ) and p(ve). 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 (α1,β1) and (α2,β2) were taken to be (1,10) and (1,10) to impose nonnegativity. The total number of samples NMCMC for the algorithm was taken to be 50,000, and the last 25,000 realizations were used to estimate statistics such as mean and standard deviation for external force f since it allows the beginning of the Markov chain. This beginning of chain is often called the burn-in. The scale parameter σX for a proposal distribution q(X*∣Xi)∝exp(-(Xi-X*)2/2σX2) can affect the efficiency of the Markov chain. In this study, σX=0.05 was used to produce a sample with reasonably low autocorrelation.Numerical results are shown in Figure5. 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 E{f} is fairly accurate compared with the exact solution, regardless of the level of noise.Numerical solutions via the proposed method.
(a)
σ
=
0.0158
(b)
σ
=
0.0316Figure6 shows examples of trace plot with respect to the number of sampling for the case σ=0.0158. These points f40 and f140 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.Examples of trace plots for the MCMC simulation and their posterior distributions.
(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 (Figure7) is generally expressed as the following nonlinear differential equation:(4.1)(I+A44)ϕ¨+B(ϕ˙)+C44ϕ=M4(t),
where I is the actual mass moment of inertia, A44 is the added mass moment of inertia, B is the nonlinear roll damping function, C44 is the coefficient of restoring moment, ϕ is the roll response, and M4 represents a wave exciting moment acting on a ship. The sum of I and A44 is often referred to as the virtual mass.Figure 7
Roll motion of a ship in waves.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.Analogy of (a) the general mechanical system and (b) the ship motion 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 mechanismC44 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 massA44 and damping mechanism B of a ship motion. The wave exciting force M4 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 excitationM4. 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. Figure9 shows the test model and its body plan illustrating the shape of hull form. Table 1 summarizes particulars of the test model.Table 1
Particulars of the model.
LengthLpp
2.500 m
BreadthB
0.387 m
DraftD
0.132 m
Displacement volume∇
0.110 m3
GM
0.074 m
Resonance frequencyωn
6.905 rad/sIllustration of (a) the test model for experimental study and (b) its body plan.
(a)
(b)Figure10 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 10
Overview of experimental setup.Figure11 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.Typical example of time history of data for a test model: wave elevation measured at far wave-probe (a) and at near wave-probe (b) and the corresponding responses (c).
(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ω=6.465rad/s in Figure 12. The data was measured with the sampling period Δt=0.01 s. Using this data, the observable quantity g in Figure 13 was calculated through (2.14). This quantity leads to the construction of the PPDF p(M4,γ,ve∣g) for wave exciting moment M4 as in (2.27).Figure 12
The measured roll response data withω=6.465rad/s.Figure 13
Observable quantityg for the response data with ω=6.465rad/s.In the same way used in the previous section, the information on unknown wave exciting momentM4 was extracted via MCMC algorithm. Without loss of generality, the values γ,ve, and M4 were initialized as γ0=1, ve0=0.1, and zero vectors for M40. The pairs of parameters (α1,β1) and (α2,β2) were also taken to be (1,10) and (1,10).Figures14 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 M4.Figure 14
Wave exciting moment obtained by pseudo inverse for the response data withω=6.465rad/s.Figure 15
Wave exciting moment obtained by stochastic inversion for the response data withω=6.465rad/s.For the purpose of validation, roll response was resimulated using the posterior mean of the identified wave exciting moment since the exact solution forM4 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 M4 for this case.Figure 16
Comparison of the re-simulated response with the measured one for the response data withω=6.465rad/s.The same procedure was also applied to the other response data whenω=3.768rad/s and ω=8.796rad/s. 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.Numerical results for the response data whenω=3.768rad/s: (a) measured response data, (b) observable value g, (c) wave exciting moment obtained by the proposed method, and (d) comparison of the re-simulated response with the measured one.
(a)
(b)
(c)
(d)Numerical results for the response data whenω=8.796rad/s: (a) measured response data, (b) observable value g, (c) wave exciting moment obtained by the proposed method, and (d) comparison of the re-simulated response with the measured one.
(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 Figure19 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 λ/L is small or large. This may be mainly because of the linear assumption in the strip theory.Figure 19
Qualitative comparison of wave exciting moment.


### 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.Figure20 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.Time history of measurements for the case of irregular wave.
(a)
(b)
(c)Figure 21
Estimated wave exciting moment through the proposed method.Figure 22
Comparison of the re-simulated response with the measured one.


## 4.1. Motion Equation of a Ship in Roll Motion

The mathematical model of ship rolling in waves (Figure7) is generally expressed as the following nonlinear differential equation:(4.1)(I+A44)ϕ¨+B(ϕ˙)+C44ϕ=M4(t),
where I is the actual mass moment of inertia, A44 is the added mass moment of inertia, B is the nonlinear roll damping function, C44 is the coefficient of restoring moment, ϕ is the roll response, and M4 represents a wave exciting moment acting on a ship. The sum of I and A44 is often referred to as the virtual mass.Figure 7
Roll motion of a ship in waves.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.Analogy of (a) the general mechanical system and (b) the ship motion 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 mechanismC44 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 massA44 and damping mechanism B of a ship motion. The wave exciting force M4 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 excitationM4. 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. Figure9 shows the test model and its body plan illustrating the shape of hull form. Table 1 summarizes particulars of the test model.Table 1
Particulars of the model.
LengthLpp
2.500 m
BreadthB
0.387 m
DraftD
0.132 m
Displacement volume∇
0.110 m3
GM
0.074 m
Resonance frequencyωn
6.905 rad/sIllustration of (a) the test model for experimental study and (b) its body plan.
(a)
(b)Figure10 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 10
Overview of experimental setup.Figure11 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.Typical example of time history of data for a test model: wave elevation measured at far wave-probe (a) and at near wave-probe (b) and the corresponding responses (c).
(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ω=6.465rad/s in Figure 12. The data was measured with the sampling period Δt=0.01 s. Using this data, the observable quantity g in Figure 13 was calculated through (2.14). This quantity leads to the construction of the PPDF p(M4,γ,ve∣g) for wave exciting moment M4 as in (2.27).Figure 12
The measured roll response data withω=6.465rad/s.Figure 13
Observable quantityg for the response data with ω=6.465rad/s.In the same way used in the previous section, the information on unknown wave exciting momentM4 was extracted via MCMC algorithm. Without loss of generality, the values γ,ve, and M4 were initialized as γ0=1, ve0=0.1, and zero vectors for M40. The pairs of parameters (α1,β1) and (α2,β2) were also taken to be (1,10) and (1,10).Figures14 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 M4.Figure 14
Wave exciting moment obtained by pseudo inverse for the response data withω=6.465rad/s.Figure 15
Wave exciting moment obtained by stochastic inversion for the response data withω=6.465rad/s.For the purpose of validation, roll response was resimulated using the posterior mean of the identified wave exciting moment since the exact solution forM4 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 M4 for this case.Figure 16
Comparison of the re-simulated response with the measured one for the response data withω=6.465rad/s.The same procedure was also applied to the other response data whenω=3.768rad/s and ω=8.796rad/s. 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.Numerical results for the response data whenω=3.768rad/s: (a) measured response data, (b) observable value g, (c) wave exciting moment obtained by the proposed method, and (d) comparison of the re-simulated response with the measured one.
(a)
(b)
(c)
(d)Numerical results for the response data whenω=8.796rad/s: (a) measured response data, (b) observable value g, (c) wave exciting moment obtained by the proposed method, and (d) comparison of the re-simulated response with the measured one.
(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 Figure19 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 λ/L is small or large. This may be mainly because of the linear assumption in the strip theory.Figure 19
Qualitative comparison of wave exciting moment.


## 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.Figure20 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.Time history of measurements for the case of irregular wave.
(a)
(b)
(c)Figure 21
Estimated wave exciting moment through the proposed method.Figure 22
Comparison of the re-simulated response with the measured one.


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



---

*Source: 175036-2012-07-18.xml*