A Matlab Toolbox for Parametric Identification of Radiation-Force Models of Ships and Offshore Structures

This article describes a Matlab toolbox for parametric identification of fluid-memory models associated with the radiation forces ships and offshore structures. Radiation forces are a key component of forceto-motion models used in simulators, motion control designs, and also for initial performance evaluation of wave-energy converters. The software described provides tools for preparing non-parmatric data and for identification with automatic model-order detection. The identification problem is considered in the frequency domain.


Introduction
One approach to develop linear time-domain models of marine structures consists of using computed hydrodynamic data for system identification and obtain a parametric approximation of the Cummins equation.If either physical-model or full-scale experiments are available, then the mathematical model based on the Cummins equation can be corrected for viscous effects.These corrections can also be obtained from computational fluid dynamics.This procedure is illustrated in Figure 1.
A great deal of work has been reported in the literature on the use of different identification methods to approximate the fluid-memory models in the Cummins equation.Taghipour et al. (2008) and Perez and Fossen (2008b) provide an up-to-date review of the different methods.In particular, the latter reference discusses the advantages of using frequency-domain methods for the identification of fluid-memory models.The data provided by hydrodynamic codes is in the frequency domain; therefore, frequency-domain identification is the natural approach to follow.This identification approach avoids transforming the data to the time domain, which, if not handled properly, can result in errors due to the finite amount of data.More importantly, frequency-domain identification allows enforcing model structure and parameter constraints.Therefore, the class of models over which the search is done is reduced, and the models obtained satisfy properties that are in agreement with the hydrodynamic modelling hypothesis.See Perez and Fossen (2008b) for further details. to perform parametric identification of fluid-memory models associated with the radiation forces of ships and offshore structures.
We address two cases.
In the first case, the infinite-frequency added mass is considered available, and the fluid-memory model is estimated based on the frequency-dependent added mass (including the infinite-frequency value) and the frequency-dependent potential damping.In the second case, the infinitefrequency added mass is considered unavailable, and it is estimated together with the fluid-memory model.That is, the complete radiation-force model is estimated.For the second case, we follow the approach proposed by Perez and Fossen (2008a).The second case is relevant for data of hydrodynamic codes based on 2D-potential theory since these codes do not normally solve the boundary-value problem associated with the infinite frequency.

Dynamics of Ships and Offshore Strucutres
The linearised equation of motion of marine structures can be formulated as The matrix M RB is the rigid-body generalised mass.The generalised-displacement vector ξ [x, y, z, φ, θ, ψ] T gives the position of the vessel with respect to an equilibrium frame (x-surge, y-sway, and zheave) and the orientation in terms of Euler angles (φroll, θ-pitch, and ψ-yaw).The generalised force vector τ [X, Y, Z, K, M, N ] T gives the respective forces in a body-fixed frame (X-surge, Y -sway, and Z-heave) and the moments about the axis of the body-fixed frame (K-roll, M -pitch, and N -yaw).This generalised-force vector can be separated into four components: (2) The first component corresponds to the radiation forces arising from the change in momentum of the fluid due to the motion of the structure and the waves generated as the result of this motion.The second component corresponds to forces due to fluid-viscous effects (skin friction and vortex shedding).The third component corresponds to restoring forces due to gravity and buoyancy.The fourth component represents the excitation pressure forces due to the incoming waves and other forces used to control the motion of the vessel.Cummins (1962) used potential theory (ideal fluid, no viscous effects) to study the radiation hydrodynamic problem in the time-domain, and found the following representation: (3) The first term in (3) represents pressure forces due the accelerations of the structure, and A ∞ is a constant positive-definite matrix called infinite-frequency added mass.The second term represents fluid-memory effects that capture the energy transfer from the motion of the structure to the radiated waves.The convolution term is known as a fluid-memory model.The kernel of the convolution term, K(t), is the matrix of retardation or memory functions (impulse responses).
By combining terms and adding the linearised restoring forces τ res = −Gξ, the Cummins Equation (Cummins, 1962) is obtained: Equation ( 4) describes the motion of ships and offshore structures in an ideal fluid provided the linearity assumption is satisfied.This model can then be embellished with non-linear components taking into account, for example, viscous effects and mooring lines-see Figure 1.

Frequency-domain Models
When the radiation forces (3) are considered in the frequency domain, they can be expressed as follows (Newman, 1977;Faltinsen, 1990): (5) The parameters A(ω) and B(ω) are the frequencydependent added mass and potential damping respectively.This representation leads to the following frequency-domain relationship between the excitation forces and the displacements: 6) Ogilvie (1964) showed the relationship between the parameters of the time-domain model (4) and frequencydomain model (6) using the Fourier Transform of (4): From expression (7) and the application of the Riemann-Lebesgue lemma, it follows that A ∞ = lim ω→∞ A(ω), and hence the name infinite-frequency added mass.
From the Fourier transform, it also follows the frequency-domain representation of the retardation functions: Expression ( 9) is key to generate data used in the identification problems that seek parametric approximations to the fluid memory term in (4).Hydrodynamic codes based on potential theory, are nowadays readily available to compute B(ω) and A(ω) for a finite set of frequencies of interest-see, for example, Beck and Reed (2001) and Bertram (2004).Hydrodynamic codes based on 3D-potential theory usually solve, the boundary-value problem associated with infinite-frequency that gives A ∞ , whereas codes based on 2D-potential theory do not normally solve this problem.

Identification of Radiation-force Models
To implement simulations models based on the Cummins equation ( 4), non-parametric fluid-memory models can be used.This method requires a discrete-time approximation of the convolution integral and saving enough past data to evaluate the convolution at each step of the simulation.This approach can be time consuming and may require significant amounts of computer memory as illustrated in Taghipour et al. (2008).
In addition, the non-parametric models are not amicable the analysis and design of vessel motion control systems.
One way to overcome these difficulties consists of approximating the fluid-memory models by a lineartime-invariant parametric model: ) where the number of components of the state vector x corresponds to the order of the approximating system and the matrices Â, B, and Ĉ are constants.Note that the above state-space approximation does not have a feed-through term D ξ in the output equation.The reason for this is that the mapping ξ → µ has relative degree 1-see Perez and Fossen (2008b) and references therein.
The approximation problem (10) can be re-casted in the frequency domain: where K(s) is matrix of rational transfer functions with entries One can estimate the transfer functions Kik (s) and then obtain the state-space model (10) via canonical realisations (Taghipour et al., 2008).The identification problem then focus on the transfer functions Kik (s).

Identification when A ∞ is Available
This problem can be formulated in terms of Least-Squares (LS) fitting: where the notation * indicates transpose complex conjugate, and w l are weighting coefficients.The nonparamtetric model K ik (jω l ) is computed via (9) using A ∞,ik and A ik (ω l ) and B ik (ω l ) for a for a given finite set of frequencies ω l .
The structure of the estimate Kik is given by ( 12), and the vector of parameters θ can be taken as From hydrodynamic properties of the model under study, it follows that the problem ( 13)-( 14) must be considered subject to the following constraints (Perez and Fossen, 2008b): The optimisation problem ( 13)-( 14) is non-linear.Two methods can be followed to solve this problem: 1. Linearise ( 13)-( 14), and solve a sequence of linear LS problems using the solution of the previous iteration to compute the weighting coefficients w l .
2. Use the solution of the linear problem to initialise a Gauss-Newton search algorithm.
The linearisation of ( 13)-( 14) is due to (Levy, 1959) and the iterative solution via a sequence of linear problems is due to (Sanathanan and Koerner, 1963): where s l,p = 1 Note that (20) results in a Linear LS minimization.After a few iterations (usually p=10 to 20), Q ik (jω l , θ p ) ≈ Q ik (jω l , θ p−1 ); and therefore, the problem ( 13)-( 14) is approximately recovered.This allows solving the nonlinear LS problem via iterations on linear LS problems.

Order Selection
With respect to the order selection of the approximation, it follows from the constraints ( 16)-( 19) that the minimum order transfer function that can be considered has the following form: Therefore, for automatic order determination, one can start with the lowest-order approximation ( 21) and increase the order to improve the fitting until a satisfactory approximation is obtained (Unneland and Perez, 2007).
As a metric for determining the quality of the fit, one can use the coefficient of determination for both added mass and damping.This coefficient is defined as where X k are the data points and Xk are the estimates.
The coefficient R 2 is interpreted as the fraction of the data that is explained by the regression model.The steps for automatic-order detection can then be as follows: (i) Obtain the parametric model ( 12).
(ii) Reconstruct the added mass and damping from the real and imaginary parts of the estimate (12), that is, (iii) Compute the corresponding coefficients of determination for added mass and damping.If either of these coefficients is below 0.99, then the model order is increased, and the procedure is repeated from step (i).Otherwise, stop.

Stability
The resulting model from the LS problem ( 13)-( 14) may not necessarily be stable because stability is not enforced as a constraint in the optimisation.This can be addressed in a sub-optimal manner.Should the obtained model be unstable, one could obtain a stable one by reflecting the unstable poles about the imaginary axis and re-computing the denominator polynomial.That is,

Passivity
The mapping ξ into a force introduced by the fluidmemory convolution is passive-see Perez and Fossen (2008b) and references therein.The frequencyresponse-LS-fitting problem ( 13)-( 14) 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.
Often, low-order approximations models given by the solution of ( 13)-( 14) are passive-the term 'low' depends on the data of the particular vessel under consideration.Therefore, one can reduce the order and trade-off fitting accuracy for passivity.

Identification when A ∞ is Unavailable
If the infinite-frequency added mass matrix A ∞ is unavailable, one cannot form K(jω) as indicated in (9).In this case, the method proposed by Perez and Fossen (2008a) can be followed and estimate jointly the infinite-frequency added mass and the fluid-memory transfer function.The method exploits the knowledge and procedures used in the identification of Kik (s) when the infinite-frequency added mass is considered available.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 On the other hand, taking the Laplace transform of (3), and assuming a rational approximation for the convolution term we obtain The transfer function in brackets in ( 29) can be further expressed as Thus, we can use LS optimisation to estimate the parameters of the approximation (30) given the frequency-respose data ( 27): with and the constraint that As already mentioned in the previous section, the minimum order approximation is n = 2. Therefore, we can start with this order and increase it to improve the fit if necessary.Hence, we can use the same algorithms that we use for the case when the infinite-frequency added mass is available, subject to different order constraints and interpretation of the estimates obtained.
If the polynomial Q ik (s) is normalised 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: 5 Toolbox Description The toolbox presented in this paper is an independent component of the Marine Systems Simulator (MSS, 2009).Figure 2 shows a diagram of the different software components of the toolbox and their dependability.The main function of the toolbox is FDIRadMod.m,which processes the input data and returns the estimate of the fluid-memory transfer function and also the infinite-frequency added mass if required.This function calls other functions to prepare the data for identification and to compute the estimates.The toolbox also includes two demos which show how to use the main function.The first demo considers the estimation with infinite-frequency added mass available (WA), and the second demo considers the estimation when infinite-frequency added mass is not available (NA).
The functionality of the main components is described in the following.

FDIRadMod.m
Purpose: This function processes the input hydrodynamic data and estimates the order and the parameters of a transfer function approximation of the fluid-memory model.The function processes only single-input-single-output models.Hence, for a multiple degree of freedom vessel or marine structure, this function should be used for each relevant coupling i, k. Syntax: Input Data: • W -Vector of frequencies.
dom (this is used to label the plots).
The structure FDIopt has the following fields: • FDIopt.OrdMax -Maximum order to be used in automatic order detection.Typical value 20.Output Data: • KrNum,KrDen -Vectors with the numerator and demoninator coefficients of the estimated single-input-single-output transfer function.• Ainf hat -Estimate of the infinite-frequency added mass coeffieicient.
If the option FDIopt.AinfFlag is set to 1, then Ainf hat=Ainf, which is part of the input data.If the option FDIopt.AinfFlag is set to 0, the input value Ainf is ignored and estimated, so the user can enter any value in the function argument.

Description: The function FDIRadMod.m first calls
EditAB.m to prepare the data for identification.Then, depending on the option FDIopt.AinfFlag, the function calls the appropriate computation routine-see Figure 2.
The function FDIRadMod.malso makes an automatic order estimate by increasing the order of the approximation and computing the coefficient of determination related to the fitting of both added mass and damping.When both these coefficients reach a value greater or equal to 0.99, the function stops increasing the order, and the reconstructued added mass and damping are plotted together with the non-parametric data used for identification.At this point, the function prompts the user to either adjust the order of the approximation manually via a keyboard input or leave the model as it is and exit the function.
The user can make as many changes in order as required, and every time there is a change in the order, the model is re-estimted and the data replotted.

EditAB.m
Purpose: This function allows the user to select the frequency range to be used for identification and to eliminate data wild points 1 .
Description: This is a support function for FDIRadMod.m,so the user may not need to call it directly.The function first plots the added mass and potential damping as a function of the frequency, and then prompts the user to select the range of frequencies for identification.This range selection is done by clicking with the mouse on the plot of either the added mass or 1 Wild points in the data computed using hydrodynamic codes are due to ill-conditioned numerical problems, which often arise at high-frequencies if inappropriate panel sizes are used to discretise the hull-see Faltinsen (1990) for details.
damping.The low-frequency data point should be selected first, then the high-frequency point, and finally the user should press the return key.
The data is then re-plotted within the selected range.
After selecting the frequency range, the function allows the elimination of data wild points.A message on the workspace prompts the user to opt for wild point elimination.If required, this elimination is done by clicking with the mouse on all the points that are to be eliminated (either on the plot of the added mass or damping), and finally the user should press the return key.The function allows the user to re-start the process in case a point is deleted accidentally.

Ident retardation FD.m
Purpose: This function estimates the parameters of a specified order approximation for the fluidmemory transfer function given the frequency response K(jω l ).
Description: This is a support function for FDIRadMod.m,so the user may not need to call it directly.This function performs the estimation for the problem in which the infinite frequency added mass is available to compute K(jω l ).This problem is described in Section 4.1.The function performs data scaling and enforces the model structure constraints ( 16)-( 19).A summary of the algorithm is given in the Appendix.

Ident retardation FDna.m
Purpose: This function performs the joint parameter estimation of the approximating fluid-memory transfer function and the infinite-frequency added mass coefficient.
Description: This is a support function for FDIRadMod.m,so the user may not need to call it directly.The function uses as input the frequency-dependant added mass and damping, and it requires a desired order.The associated estimation problem is described in Section 4.5.
The function performs data scaling and enforces the model structure constraints ( 16)-( 19).

Fit siso fresp.m
Purpose: This is a general purpose function to estimate a single-input-single-output transfer function of a specified order and relative degree given a frequency response.
This function can be used not only to identify fluid-memory transfer functions, but also forceto-motion transfer functions (Perez and Lande, 2006).That latter is a functionality that will be included in future versions of the toolbox.
Description: This is a support function for Ident retardation FD.m and Ident retardation FDna.m, so the user may not need to call it directly.
This function implements 3 methods for parameter estimation, namely, 1 -uses a linearised model and linear LS optimisation.2 -uses iterative linear LS optimisation, 3-uses the linear LS optimisation solution to initialise a non-linear LS optimisation problem solved using the Gauss-Newton method.The function is built upon the functionality invfreqs.m of Matlab's Signal Processing Toolbox.

Demos
The toolbox provides two demo files that make use of the main function FDIRadMod.m-seeFigure 2.These demos are based on the data of a FPSO that belongs to the Hydro add in of the Marine Systems Simulator (MSS, 2009).
The first demo, Demo FDIRadMod WA.m (WA-with infinite-frequency Added mass), loads the data structure vessel, and allows the user to select the desired coupling (i, k) for identification.The structure vessel contains data corresponding to 6 degrees of freedom, that is, i, k = 1,...,6.In this section, we illustrate the estimation results on the models corresponding to vertical motion modes; that is, couplings 3-3, 3-5, 5-3, and 5-5.
Figure 3 shows the raw added mass and damping for coupling 5-3, which by symmetry of the hull it is the same as the 3-5 coupling.These data are obtained from a hydrodynamic code.Figure 4 shows the edited data after eliminating some wild points.Figure 5 shows the corresponding curve fitting results.This figure shows the fitting of the fluid-memory frequency response on the left-hand side and the re-construction of added mass and damping based on the estimated model on the right-hand side.The order of the approximation is 5, which is obtained automatically by the function.Figures 6 and 7 show the corresponding results for the 3-3 and 5-5 couplings.For both these couplings the automatic order detection selected order 3, but then we manually increase the order to 4 to improve the fit.
The second demo, Demo FDIRadMod NA.m (NA-No infinite-frequency Added mass), also loads the data structure vessel of the FPSO, and allows the user to select the desired coupling for identification.In this demo, however, the identification is done without using the infinite-frequency added mass coefficient.Figure 8 shows the fitting results for the coupling 5-3.The left-hand-side plots show the fitting of the complex coefficient Ã(jω) given by ( 27), whereas the righthand-side plots show the re-construction of added mass and damping based on the estimated model. Figure 9 shows the estimated fluid-memory frequency response function.These results are in agreement with those shown in Figure 5; however, there are small differences due to the fact that the two estimators use different information.
It is worthwhile highlighting that the removal of wild points could be very important.For example, Figures 10 and 11 show the results of identification without using added mass for the coupling 5-3 when the wild points in the added mass and damping have not been removed.In this case, the automatic order detection selected an approximation of order 10.This is because the algorithm tries to fit lightly-damped complex poles to the wild points.The function gives the user the option to manually reduce the order.However, for some cases this may not solve the problem, and the identification process should be started again and remove the wild points.

Dependence on Other Matlab Toolboxes
The toolbox presented in this paper is based on standard Matlab code, except for two specific functions that belong to the Signal Processing Toolbox, namely, freqs.m and invfreqs.m.The function freqs.m computes the frequency response of a transfer function model for a specified set of frequencies.The function invfreqs.mperforms the parameter estimation using either the linear LS method or the Gauss-Newton method.In order to implement the iterative linear LS method, this function is called recursively by fit siso fresp.m.
If the user does not have a licence for the Signal Processing Toolbox, the two functions used should be coded by the user.An alternative to the function freqs.m is very simple to code.The algorithm given in the Appendix can be used as a guide to implement an alternative to invfreqs.m.

Software Repository
The toolbox presented in this paper is an independent component of the Marine Systems Simulator (MSS, 2009) maintained by the authors.This is a free toolbox released under a GNU licence.The software is available at www.marinecontrol.org

Conclusion
This paper describes a toolbox for parametric identification of fluid-memory models associated with the radiation forces of marine structures based on a frequency-domain method.The models identified find application in the development of ship simulators, control design, and the evaluation of wave energy converters.
The software described provides tools for preparing the non-paramatric data generated hydrodynamic codes, automatic model order detection, and parameter estimation.The toolbox contains a main function that performs all these tasks by calling other support functions.The user may only need to call the main function.
The identification is done for single-input-singleoutput models.This gives the user freedom to select the couplings of interest each particular application and to integrate the functionality of the toolbox into other data processing codes.

Figure 3 :Figure 4 :
Figure 3: Raw added mass and damping data of a FPSO vessel computed by a hydrodynamic code.Coupling 5-3 (pitch-heave).

Figure 5 :
Figure 5: Identification results for the coupling 5-3 (pitch-heave) using information of the infinite-frequency added mass.The left-hand-side plots show the fluid-memory frequency response data and the response of the identified model.The right-hand-side plots show the added mass and potential damping and their re-construction from the estimated model.

Figure 6 :
Figure 6: Identification results for the coupling 3-3 (heave-heave) using information of the infinite-frequency added mass.The left-hand-side plots show the fluid-memory frequency response data and the response of the identified model.The right-hand-side plots show the added mass and potential damping and their re-construction based the estimated model.

Figure 7 :
Figure 7: Identification results for the coupling 5-5 (pitch-pitch) using information of the infinite-frequency added mass.The left-hand-side plots show the fluid-memory frequency response data and the response of the identified model.The right-hand-side plots show the added mass and potential damping and their re-construction based on the estimated model.

Figure 8 :
Figure 8: Identification results for the coupling 5-3 (pitch-heave) without using information of the infinitefrequency added mass.The left-hand-side plots show the complex coefficient Ã(jω) data and the response of the identified model.The right-hand-side plots show the added mass and potential damping and their re-construction based on the estimated model.

Figure 9 :
Figure 9: Frequency response of the identified fluid-memory model for the coupling 5-3 (pitch-heave) without using information of the infinite-frequency added mass.

Figure 10 :
Figure 10: Identification results for the coupling 5-3 (pitch-heave) without using information of the infinitefrequency added mass and without eliminating wild points.The left-hand-side plots show the complex coefficient Ã(jω) data and the response of the identified model.The right-hand-side plots show the added mass and potential damping and their re-construction based on the estimated model.

Figure 11 :
Figure11: Frequency response of the identified fluid-memory model for the coupling 5-3 (pitch-heave) without using information of the infinite-frequency added mass and without eliminating wild points.

•
FDIopt.AinfFlag -Logic flag.If set to 1, the value Ainf is used in the calculations.If set to 0, the infinite-frequency added mass is estimated, and the value in the argument of the function is ignored.