Joint Identification of Infinite-Frequency Added Mass and Fluid-memory Models of Marine Structures

This paper addresses the problem of joint identification of infinite-frequency added mass and fluid memory models of marine structures from finite frequency data. This problem is relevant for cases where the code used to compute the hydrodynamic coefficients of the marine structure does not give the infinite-frequency added mass. This case is typical of codes based on 2D-potential theory since most 3D-potential-theory codes solve the boundary value associated with the infinite frequency. The method proposed in this paper presents a simpler alternative approach to other methods previously presented in the literature. The advantage of the proposed method is that the same identification procedure can be used to identify the fluid-memory models with or without having access to the infinite-frequency added mass coefficient. Therefore, it provides an extension that puts the two identification problems into the same framework. The method also exploits the constraints related to relative degree and low-frequency asymptotic values of the hydrodynamic coefficients derived from the physics of the problem, which are used as prior information to refine the obtained models.


Introduction
Time-domain models for rigid-body motion simulation of marine structures are of paramount importance for the development of training simulators, hardware-inthe loop testing simulators, wave energy converters and motion control systems.One way to develop these models consist of using hydrodynamic codes based on potential theory to compute frequency frequencydependent coefficients, and then use these data to obtain time-domain models via system identification.Two approaches can be followed for the latter part.One approach consists of using the Cummins equation (Cummins, 1962).The other approach consists of using the force-to-motion data directly (Perez and Lande, 2006).In this paper, we concentrate on the first approach.
The Cummins Equation relates the motion of the marine structure to the wave-induced forces in time domain under the assumption of linearity.This equation is an integro-differential equation that contains a convolution term representing fluid-memory effects associated with the dynamics of the radiation forces.This convolution term is inconvenient for simulation and also for the analysis and design of motion control systems.A linear-time-invariant model can be used to approximate the convolution in the Cummins equation.To obtain such a linear model, one can apply system identification.The identification problem can be posed either in the time or in the frequency domain.Due to these alternative problem formulations, there has been a great deal of work reported in the literature-see, for example, Jefferys et al. (1984), Jefferys and Goheen (1992), Yu and Falnes (1995), Holappa and Falzarano (1999), Hjulstad et al. (2004), Kristansen and Egeland (2003), Kristiansen et al. (2005), Jordan and Beltran-Aguedo (2004), McCabe et al. (2005), and Sutulo and Guedes-Soares (2005).Taghipour et al. (2008) provide a review and a summary of some of the methods.Perez and Fossen (2008) compared and discussed the advantages and disadvantages of time-and frequency domain methods for the identification of the fluid memory models.It is argued that the frequency-domain identification results in estimation algorithms that are easier to implement and use than those resulting from time-domain formulations.Also, the quality of the models obtained is, in general, superior to those obtained with the time-domain methods proposed in the literature.The simplicity of the identification stems from the fact that the solution to the parameter estimation problem can be based on iterative linear Least Square (LS) optimisation.Prior knowledge derived from the hydrodynamics of the problem can be used to set constraints on the parameters, and the use of constraints in LS estimation leads, in general, to more accurate estimates-see e.g., (Gourieroux and Monfort, 1995).On the negative side, It has also been discussed that the frequency-domain identification approach can be sensitive to the estimate of the infinitefrequency added mass coefficient provided by the 3Dhydrodynamic codes.The need of a reliable estimate of the infinite-frequency coefficients poses a bigger problem when using 2D-hydrodynamic codes, since these codes do not provide such estimate.
In this paper, we present a procedure that extends the application of frequency-domain identification of seakeeping models to the case where the hydrodynamic data does not include the infinite-frequency added mass coefficient (or one choses not to use it).That is, only finite frequency data is considered.The model identified relates the total radiation forces to the velocities, and it allows identifying the infinite-frequency added mass together with a fluid memory model.The proposed method is motivated by the work of Kaasen and Mo (2004), but provides a simpler alternative.

A Linear Model based on Cummins Equation
The equations of motion of a rigid marine structure in body-fixed coordinates can be linearised about an equilibrium point and heading and be expressed as where ξ represents the generalised perturbation position-orientation vector, τ is the vector of generalised forces and moments, and M RB is the positivedefinite rigid-body generalised inertia matrix.The generalised force vector τ can be separated into three components: where the first component corresponds to the radiation forces arising from the change in momentum of the fluid due to the motion of the structure, the second are restoring forces due to gravity and buoyancy, and the third component represents the pressure forces due to the incoming waves.For further background information about these models see for example Newman (1977) and Faltinsen (1990).Cummins (1962) studied the radiation hydrodynamic problem in an ideal fluid in the time-domain and found the following representation for the linear pressure forces: (3) The first term in (3) represents forces due the accelerations of the structure, and A is the constant positive definite added inertia matrix.The second term represents fluid memory effects that incorporate the energy dissipation due the radiated waves consequence of the motion of the structure.The kernel of the convolution term, K(t), is the matrix of retardation or memory functions (impulse responses).By renaming the variables and combining (1), (2), and (3), we obtain the Cummins Equation: Equation (4) describes the motion of the structure for any wave excitation τ exc (t) provided the linearity assumption is satisfied; and it forms the basis of more complex models, which can be obtained by adding nonlinear terms to represent different physical effects.
There are zeros at s = 0.

Frequency-domain Representation of the Radiation Forces
When (3) is considered in the frequency domain, it takes the following form (Newman, 1977;Faltinsen, 1990): The parameters A(ω) and B(ω) are the frequencydependent added mass and damping respectively.Hydrodynamic codes based on potential theories (2D and 3D) are nowadays readily available for the computation of the frequency-dependent added mass A(ω) and potential damping B(ω).These data are computed for a reduced set of frequencies of interest.Ogilvie (1964) showed using the Fourier Transform of (4), that the following frequency-domain representation holds for the retardation functions: and also that from which the name infinite-frequency added mass follows.

Frequency-domain Identification of the Convolution Terms
Expression (6) provides a way to compute the frequency response function K(jω) for a finite set of frequencies.These data is the basis for the frequencydomain identification methods that seek a transfer function approximation to each entry of K(jω): Apart from the non-parametric frequency-response data K ik (jω), there is prior information that should be used as much as possible to refine the search for the appropriate model and its parameters.This is an important aspect of any identification problem since, in general, using prior information to set constraints on the model structure and parameters leads to more accurate estimates.Table 1 summarizes the properties of the retardation functions and their implications on the parametric models (8).The left column shows properties in frequency-and time-domain that derive from the hydrodynamics of the problem under considerations.For example, the first and second property are the are consequence of no waves being generated due to the motion of the structure at zero and infinite frequency.The third property derives from the fact that K(0 + ) equals the area under the curve of B(ω).The fifth property is related to the dissipative characteristics of the radiation forces.For further discussion of these properties and their derivations see Perez and Fossen (2008) and references therein.
The properties shown on the left column of Table 1 have consequences on the models (8), and these are shown in the right column of the table.These properties are related to the structure of the models.Indeed, we can express the models (8) as From Table 1, it is known that these transfer functions have a zero at s = 0, hence, this information can be taken into account by writing the models as with the constraint on the order of the polynomials Since, in general, A(0) = A(∞), it follows from (6) that there is unique zero of K ik (s) ar s=0.Therefore, l = 1, and m = n − 2. This is simple to verify from the non-parametric data since the phase of K ik (jω) at low frequencies tends to l π/2.
One way to exploit this information in the identification process is to consider as data for the identification of P ′ ik (s) and Q ik (s) with the constraints The identification problem can be posed as a complex curve fitting problem: with and the vector of parameters , θ is defined as θ = [p m , ..., p 0 , q n−1 , ..., q 0 ] T .
The weights w l can be exploited to select how important is the fit at different frequency ranges.The above parameter estimation problem is a nonlinear LS problem in the parameters, which can be solved using a Gauss-Newton algorithm, or it can be linearized as indicated in the next section.Levy (1959), proposed a linearisation of ( 13)

A Linear Iterative Solution
with This problem is linear in the parameters, and thus easy to solve.Indeed, using a matrix form we can write with Using this notation, we can write with the obvious definition for the matrices Φ and Γ.
The solution to ( 18)-( 20) is then given by The linearised problem ( 16) derives from the non-linear problem (13) by choosing This means that solving ( 16) can be thought of as solving ( 13) with the weights as given in ( 22).The weights w ′ l normally correspond to a rectangular window.
A problem with this linear formulation is that the identified transfer function does not in general give a good fitting.For example, when the data extends over a large range of frequencies or when Q(s) has poorly damped complex roots close to the imaginary axis, the weighting coefficients w l in (22) will weight the fit more heavily at high frequencies, and also at frequencies close to the resonant roots.This may give a bias in the parameter estimates.Sanathanan and Koerner (1963) proposed a method to compensate for the bias introduced by the linearisation.This method consists in solving ( 18)-( 20) iteratively using as weighting coefficients the inverse of the denominator Q(jω, θ) evaluated at the previous estimate.This algorithm can be summarised in the following: 1. Set W 0 = I.

Solve θ
This choice of weighting coefficients in step 3 results in the following problem at each step k of the iteration: Normally, after a few iterations (10 to 20), θ ⋆ k ≈ θ ⋆ k−1 ; and thus, the original non-linear LS problem ( 13) is approximately recovered.

Oder selection
The order of the transfer functions depends on the hydrodynamic characteristics of the vessel; i.e., it depends on the hull shape.Based on the properties of the convolution terms given in Table 1, it follows that the minimum order transfer function that satisfies all the properties is a second order one: Therefore, one can start with this minimum order transfer function (n=2), and increase the order while monitoring that the LS cost decreases-or simply by visual inspection of the fitted frequency response.If the order of the proposed model is too large, there will be over-fitting and therefore, the cost will increase; however before this happens, the value of the cost normally remains unchanged as one increments the order of the system.

Stability
The resulting model from the LS minimization may not necessarily be stable because stability is not enforced as a constraint in the optimisation.This can be addressed by reflecting the unstable poles about the imaginary axis and re-computing the denominator polynomial.That is, • Compute the roots of λ 1 , . . ., λ n of Q ik (s, θik ).

Passivity
It also follows from the properties given in Table 1, that the diagonal terms K ii (jω) are passive; i.e., the real part B ii (ω) must me positive for all frequencies.
For the off-diagonal terms K ik (jω) this may not be the case however.The method of LS curve fitting is that it does not enforce passivity.If passivity is required (i.e., B ik (ω) > 0), a simple way to ensure it is to try different order approximations and choose the one that is passive.The approximation is passive if When this is checked, one should evaluate the transfer function at low and high frequencies-below and above the frequencies used for the parameter estimation.
Normally, the low-order approximations models of the convolution terms given by this method are passive.Therefore, one can reduce the order and trade-off fitting accuracy for passivity.A different approach would be optimise the numerator of the obtained non passive model to obtain a passive approximation-this goes beyond the scope of this paper, but the reader is referred to Damaren (2000) and references therein.
3 Joint Identification of Infinite-frequency Added Mass and Fluid Memory Models Hydrodynamic codes based on 2-D potential theory normally do not provide the value of the infinite frequency added mass coefficient A = lim ω→∞ A(ω).In these cases, we cannot form K(jω) as indicated in ( 6).Kaasen and Mo (2004) addressed this problem by making a partial-fraction expansion of the real part of Kik (s) in terms of ω 2 .The poles and residuals of this expansion can be estimated using Least-squares and the damping data B ik (ω).Then Kik (s) can be obtained by mapping poles and the residuals of the partial-fraction expansion of its real part into the poles and residuals of a partial fraction expansion of Kik (s) in terms jω.
In this section, we propose a simpler alternative to the method of Kaasen and Mo (2004).The proposed method exploits the knowledge and methods used in the identification of Kik (s) discussed in the Section 2.2, and therefore, it provides an extension of those results putting the two identification problems into the same framework.
On the one hand, the radiation forces in the frequency-domain given in (5) can be expressed where the expression in brackets gives the complex coefficient Ã(jω) On the other hand, taking the Laplace transform of (3), and assuming a rational approximation for the convolution term we obtain This representation can be traced back to the work of Söding (1982), and it has been used by Xia et al. (1998) and Sutulo and Guedes-Soares (2005), but with a different approach to that presented in this paper.The transfer function in brackets in (28) can be further expressed as Thus, we can follow the same approach as in Section 2.2 and use Least-squares optimisation to estimate the parameters of the approximation (29) given the frequency-respose data (26): with and the constraint that n = degR ik (s) = degR ik (s).
We also know from Section 2.4 that the minimum order approximation is of n = 2. Therefore, we can start with this order and increment to improve the fit if necessary.
It should be noted as well that since we have normalised the polynomial Q ik (s) to be monic, then that is, the infinite-frequency added mass A ik is the coefficient of the highest order term of R ik (s, θ ⋆ ).Also, after obtaining R ik (s, θ ⋆ ) and S ik (s, θ ⋆ ), we can recover the polynomials for the fluid-memory model: 4 Model Quality Assessment In order to assess the quality of the model, one can compare the frequency-dependant added mass A ik (ω) and damping coefficients B ik (ω) provided by the hydrodynamic code, with those reconstructed from the estimated retardation function: Good fitting of these coefficients give confidence in the estimated values of Kik (s) and Âik .

Case Studies
To illustrate the use of the method proposed in the previous section, we consider the hydrodynamic data of three different vessels: • Containership, • FPSO, • Semi-submersible.
The hydrodynamic data for all vessels is computed with WAMIT.This code gives an estimate of the value of the infinite-frequency added mass; and therefore, we have means to validate the estimated parameters if we assume the values given by the code are close to the true parameters.The containership vessel is the same vessel used in (Taghipour et al., 2008) (Perez and Fossen, 2008).The FPSO and the Semi-submersible are the example demos provided with the Marine Systems Simulator (www.marinecontrol.org).
For the containership vessel, we have computed only the part of the model corresponding the vertical-plane motion (heave and pitch).Table 2 shows the results of the estimated infinite-frequency added mass coefficients, together with the true values and the absolute relative error.Figure 1 shows the fitting of the complex coefficient Ã33 (jω), the reconstruction of the added mass and damping based on (34).As we can see from this figure, the fitting is relatively good.
Tables 3 and 4 show the estimation results in six degrees of freedom for the FPSO and Semi-submersible respectively.Note that since the semi-submesible hull has fore-aft and port-starboard symmetry, there are less couplings.Figures 2, 3, and 4 show the fit for particular couplings.
As we can see from the examples used in this section, the method is able to estimate the infinite-frequency added mass coefficient with good accuracy and also to provide high-order fittings as those shown for semisubmersible.

Conclusions
This paper addresses the problem of joint identification of infinite-frequency added mass and fluid memory models from finite-frequency data.This problem is particularly relevant to the cases where the hydrodynamic code used to compute the coefficients does not give the infinite frequency added mass coefficient.This is the case for codes based on 2D potential theory.This problem has been previously addressed via partial-fraction expansions by Kaasen and Mo (2004).The method proposed in this paper presents simpler A 11 = 7.363e06 Â11 = 7.546e06 2.5 % A 22 = 3.393e07 Â22 = 3.589e07 5.8 % A 33 = 5.929e07 Â33 = 6.023e07 1.6 % A 44 = 6.065e10Â44 = 6.083e10 0.3 % A 55 = 5.021e10 Â55 = 4.975e10 0.9 % A 66 = 3.756e10 Â66 = 3.729e10 0.7 % A 13 = 8.663e07 Â13 = 9.570e7 10 % A 24 = -4.624e08Â24 = -5.075e89.7 % alternative to the existing proposal.The advantage of the method is that the same identification procedure can be used to identify Kik (s) when the A ik is given and Kik (s) and Âik when the latter is not given.The method also exploits the information related to relative degree and low-frequency asymptotic values of the hydrodynamic coefficients derived from the physics of the problem.This information is used to impose constraints on the model structure.

Figure 1 :
Figure 1: Fitting results for the containership.Left column: Frequency response Ã33 (jω) and estimate.Right column: Reconstruction of added mass and damping from the identified fluid memory function K33 (jω) based on an 3rd order approximation.

Figure 2 :
Figure 2: Fitting results for the FPSO.Left column: Frequency response Ã33 (jω) and estimate.Right column: Reconstruction of added mass and damping from the identified fluid memory function K33 (jω) based on an 2nd order approximation.

Table 1 :
Properties of Retardation Functions PropertyImplication on Parametric Models K ik

Table 2 :
True and Identified Added Mass Coefficients