A Ship-handling Simulator
Sea Wave Modeling 

SADKO MANDŽUKA, ACO ŠIKANIĆ, DARIO SEFEROVIĆ
Brodarski Institute, d.o.o.
Av.V.Holjevca 20, 10000 Zagreb
CROATIA

Abstract: - The different mathematical models are implemented for ship handling simulators.  The main parts of ship handling simulator are: mathematical model of own ship (library of some ship models), mathematical model of environment (waves, wind and current), various port layouts etc. The characteristics of sea disturbances have to be defined quantitatively during the design of ship simulators (ship handling, ship navigation, marine engine simulators, etc.).  The convenient data are available in the form of the energy spectral density functions (Pierson-Moskowitz, Jonswap, Voznesensky-Netsvetaeva, Darbyshire, Tabain, Davidon). Following a linear theory, using the spectral factorization techniques, the sea spectrum has been modeled by transfer function representation (colored filter). An algorithm of non-linear regression for rationalization of sea spectrum, described by the Pierson-Moskowitz expression, is given in the paper. The results of the computer simulation are presented.

Key-Words:  - ship-handling, sea wave modeling, sea spectra, colored filter

1   Introduction

	There are at least two cases when the modeling of sea spectra is needed. Firstly, it is necessary for digital simulation of the behavior of floating vessels. To analyze more complex systems on floating vessels and verify the design usually requires a simulation. These simulations use a computer to solve some algebraic and differential equations which model the system. The wave profile was at first generated by approximation of the wave spectrum with a finite number of discrete frequency sinusoids with random phases, Morgan, [1]. The number of sinusoids used to approximate the wave spectrum varies depending on the desired accuracy in simulating the spectral density function selected to model the real-world waves. Using too many sinusoids to simulate the wave spectrum is of questionable value. The method has a pronounced disadvantage if the simulation is being performed with a computer of limited capabilities. It is natural to move away from this approach towards modern stochastic processes theory which is based on a more accurate and simplistic methodology.
	Secondly, the sea disturbances have to be defined adequately in the process of control systems design using modern optimal control and estimation theory, Papoulis [2], Grimble [3]. The traditional approach, using the classical seakeeping methodology, Bishop and Price [4] and Lloyd [5], is not suitable in this 

case. Lately, the progress obtained in the stochastic process theory has been very remarkable, not only in mathematics and physics, but also in communications, measuring and control engineering, as well as economics and biology. The modern control and estimation theory are mostly based on linear models. The assumption which is fundamental to the development of behavioral models of the vessels on waves for control design is that the sea state is known and can be described by a spectral density function. It is assumed that in the worst case, the high frequency motions of the vessel are the same as the oscillatory motions of the waves, Mandžuka [6]. Consequently, the behavior of the vessel in response to these high frequency forces may be modeled by colored filters for each degree of freedom of motion.
	The concept of the colored filter realization of sea spectrum has been adopted rather early in the control and estimation theory for ships. Unfortunately, the incorrect results are given in some papers, Li and Leigh [7], in some studies inadequate models are used, and in most of them the mathematical foundation for this approach is not described, Van Amerongen, [8].
	The methodology of color filter synthesis is described in Section 2. The initial assumption for that synthesis is that the sea spectrum is given in a rational proper function form of square frequency (Markov process). This means that the roots of numerator and denominator are symmetric with respect to the real and imaginary axes. The transfer function of color filter can then be chosen so that all its poles are in the left half-plane and all its zeros in the left half-plane or on the imaginary axis. This technique is known as spectral factorization, Демиденко, [9].
	The first step in color filter design is in spectrum approximation in the form of rational proper function of square frequency. This is the classical non-linear least squares problem, or general data fitting, Fletcher [10]. In Section 3, a new optimization algorithm is suggested. It is based on transformation of originally non-linear problem to a recursive linear form.
The optimization algorithm is realized in the MATLAB (Mathworks, Inc). In Section 4, the values for the coefficients required to set up the filters are tabulated as a function of significant wave heights. The graph presentations of wave spectrum and their approximations are given. It is shown that the agreement is very good in the large interval around resonant frequency.

2    Color Filter Synthesis

A link of an input and output spectrum of a dynamic system is a well-known result in stochastic processes theory, Papoulis [2], and can be given by:

	(1)	

If the excitation signal x(t) is white noise with characteristics: E{x(t)}=0; E{x2(t)}=1 (unit white noise), and if the output spectrum Syy(ω) is given in a rational proper transfer form of square frequency, then the transfer function F(s), also called color filter (shaping filter), can be realized. The output spectrum can then be given in the form:

	(2)	
where:
	(3)	
If pi and qi are real numbers, then the roots of corresponding polynomials are either real or complex conjugate. The P(ω) and Q(ω) polynomials can then be written in the form:
	(4)	
where Qpk and Qqk are all the roots above the real axis and pm1 and qm2 are the polynomials' coefficients of the highest power. The conjugated values are marked by (*). If each of polynomials (4) is repeatedly multiplied by (j)*(-j) =1, where j is imaginary unit, then these polynomials can be written in the form:

	(5)	

The singularities of complex spectrum are symmetrical with respect to the jω axis and can be separated into two groups. The "left" group consists of all si with Re{si} 0. This is causal part. The "right" group consists of all other singularities. Then, the polynomials (5) can be written in symbolic form:
	(6)	

where all the roots of p+(jw) and q+(jw) lie in the left half-plane (minimum phase component) . Since:
	(7)	
the corresponding transfer function is:
	(8)	
or :
	(9)	
where s is Laplace operator.
Apart from this problem, an inversion problem may also be given. It is the filter synthesis problem of white noise generating if the signal of known spectrum density power is brought to input. Its transfer function is equal to a reciprocal transfer function of the colored filter.
In this way, the procedure for getting the colored filter if the spectrum is given in rational proper transfer form of square frequency is obtained.

3. Spectrum Rationalization

The first step in colored filter design is in spectrum approximation in the form of rational proper function (of square frequency). The problem may be given in the form of minimizing the integral squared error cost function:
	(10)	

where ω1 and ω2 define frequency interval (window) of interest, and F(ω) is the wanted approximation in the form:

	(11)	

Set b0=1 to have unique solution. The abbreviated form of approximation (11) can be written as:
	(12)	
where:

If the frequency window is discrete at n points, then the optimization criterion can be written in discrete form:
	(13)	
where:

	(14)	

and ρi is a weighting factor that defines the window shape, that is, the frequencies of interest can be weighted. The problem presented in this way, represents the classical non-linear least squares problem, or data fitting. The application of well known conventional numerical procedures (Gauss-Newton method or generalized least squares method, Nelder-Mead simplex method, etc.) has highlighted various problems such as slow convergence, rank deficiency and similar. Due to that reason, the modified cost function in recursive form is 

suggested:
	(15)	

for which it is valid:

	(16)	
By rearranging the recursive cost function, it is finally obtained:
	(17)	

By this mathematical trick the original non-linear problem is transformed into linear least squares problem. Defining the cost function in that "ad hoc" way, the analytical solution in k-th iteration is made possible. Namely, the following is obtained from the optimum conditions:

	(18)	

After rearranging the above equations, the following is obtained:

	(19)	

The system of equations (19) can be written in the matrix form:
	(20)	

with the separate partial sums in submatrix form:

	(21)  	

The matrices X and Y are formed from vectors xi and yi as their columns (dimensions are m1*n and m2*n, respectively), the matrices P and S are diagonal formed of elements ρi* and Si, and  the elements of vector v with dimension n*1 are all ones. By solving the system (20) the solution in k-th iteration can be obtained. The procedure is repeated with a change b[k+1]=b[k] until the accuracy is achieved. The initial value for b[0] is recommended to be the unit vector. The important feature of any recursive algorithm is the convergence test which is required to terminate the iteration (stopping rule). In this case, it would be demanded that, where ε is user-supplied parameter.
4.  Optimizations results
	The first step in identification of colored filter transfer function is in determination of its structure. Consequently, there exists a classic compromise problem between approximation quality (accuracy) and simplicity (complexity) of a filter. Common result of the majority of research efforts in the area of sea spectrum modeling can be summarized with the following structure:
	(22)	
It represents agreeable compromise solution for real 

range of sea states. In this case, the maximum slope in low-frequency range is defined by the numerator (40 dB/dec), while the maximum slope in high-frequency range is determined by difference of degrees of numerator and denominator (-40 dB/dec). More complex spectrum shapes indicate the more complex transfer function structure. 
	The optimization algorithm is realized in MATLAB. There are two parts of computer program: the approximation of Pierson-Moskowitz spectrum by rational proper function and colored filter design. The input data are: significant wave height, frequency limits (or dynamic range, optionally), weighting factor (default value is one). In the case that an unstable solution for F(s) is obtained, a correcting routine is used to reflect the offending poles around the imaginary axis. Very good convergence is demonstrated, and the solution for stopping rule:
ε  = 10-8
is obtained already after 5-6 iterations. The results for a few significant wave heights are shown in Table 1. Sometimes it is desirable that the dispersion of the output of colored filter is equal to the dispersion of original spectrum. 

hs           	0	1	2	a1	a2	 a3 	a4	K	
1.1	1.197	0.661	9.700	1.442	4.568	2.435	3.474	0.235	
 2.0	0.888	0.490	7.194	1.069	2.512	0.993	1.051	0.273	
3.1	0.713	0.394	5.778	0.859	1.621	0.515	0.437	0.305	
4.5         	0.592	0.327	4.796	0.713	1.117	0.294	0.208	0.334	
  6.7	0.485	0.268	3.930	0.584	0.750	0.162	0.094	0.369	
  9.2	0.414	0.229	3.354	0.499	0.546	0.101	0.050	0.400	
 12.3	0.358	0.198	2.901	0.431	0.409	0.065	0.028	0.430	
15.5	0.319	0.176	2.584	0.384	0.324	0.046	0.018	0.455	
	
Table 1. Colored filter coefficients

The graph presentation of wave spectrum (solid line) and its approximation (dotted line) is given in Figure 1. It is obvious that the agreement is very good in the greatest interval around resonant frequency.
5.   Conclusion
	The presented procedure of wave spectrum rationalization significantly improves the possibilities of application of linear mathematical models for floating vessels behavior. Using the recurrent form of the cost function many problems which are usually present in non-linear fitting can be avoided. Based on these results, it is possible to make the corresponding tables of the transfer function coefficients depending on significant wave height.
	The future research will deal with the determination of function dependence of some transfer function coefficients related to significant wave height (regression curve). The relations between statistical characteristics of sea state (significant wave height) and transfer function coefficients could be obtained. Apart from this application, the presented procedure may be used in any other case when the energy spectrum is known and anticipated to be presented in the form of the corresponding colored filter (acoustical spectrum of a ship, wind spectrum, etc.).

Fig 1. Original spectrum and its approximation

Acknowledgements

	This work is supported by a grant from the Ministry of Science, Education and Sport of the Republic of Croatia (Contract 942-20-11-2003). The presented results are developed in the framework of the Information Technology Project 2003-057.

References:

[1 ]	Morgan, M., 1978, Dynamic Positioning of Offshore Vessels, Division of The Petroleum Publishing Company, Tulsa, Oklahoma, USA.
[2 ]	Papoulis, A., Probability, Random variables and Stochastic Processes, McGraw Hill, NewYork, 1965, 342-350. 
[3 ]	Grimble, M.J., Patton, R.J., Wise, D.A., Use of Kalman filtering techniques in dynamic ship-positioning systems, IEE Proceedings, Vol. 127., No. 1, 1980, 93-102.
[4 ]	Bishop, R.E.D., Price, W.G., 1974, Probabilistic theory of ship dynamics, Chapman & Hall, London. 
[5 ]	Lloyd, A.R.J.M., SEAKEEPING: Ship Behaviour in Rough Weather, JW&S, Chichester, 1989, 106-112. 
[6 ]	Mandžuka, S., Some characteristics of sea spectrum modelling by coloured filter, 833-841, Proceedings of International Symposium: Waves - Physical and Numerical Modelling, Vancouver, 1994.
[7 ]	Li, C.K., Leigh, J.R., 1983, S-domain realisation of sea spectrum, Electronics letters, Vol. 19, No. 5, 186-187.
[8 ]	Van Amerongen, J., 1984, Estimating the wave spectrum and the roll angle and its derivatives with an extended Kalman filter, Proc. 7th Ship Control System Symp., Vol. 4, 73-85, Bath.
[9 ]	Демиденко, Е.З.,Оптимизация и регрессия, Наука, Москва, 1989, 216-232. 
[10 ]	 Fletcher, R., 1987, Practical methods of optimisation, John Wiley and Sons, Chichester.
[1 ]