J Mar Sci Technol (2007) 12:22-33 D O I 10.1007/s00773-006-p239-0
O R I G I N A L ARTICLE
Simulation of the mean zero-up-crossing wave period using artificial
neural networks trained with a simulated annealing algorithm
Haniid Bazaigan • Hainid Bahai • Faizad Aryana
Received: August 9, 2005 / Accepted: November 7, 2006 © JASNAOE 2007
Abstract The aim of this work was to develop a predic-tive model to forecast the mean zero-up-crossing wave periods (T,) for 3-hourly sea states at a location i n the Pacific using artificial neural networks (ANNs). Seven multilayer A N N s were trained with a simulated anneal-ing algorithm. The output of each trained A N N was used to estimate each of the seven parameters of a new distribution called the hepta-parameter spline proposed f o r the conditional distribution of T^, given some mean zero-up-crossing wave periods and significant wave heights. After estimating the parameters of the distribu-tion, the model was used to simulate and predict future values of 21. Forecasting a sea state and developing the joint distribution of sea state characteristics with the
help of the simulated characteristics are also discussed in this article.
Key words Mean zero-up-crossing wave period • Simulation • Neural networks • Simulated annealiitg • Hepta-parameter spline distribution
List of symbols
Cll a parameter of the hepta-parameter spline distribution
H . Bazargan • F. Aryana
College of Engineering, Shahid Bahonar University of Kerman, Iran
H . Bahai ( H )
School of Engineering and Design, Brunei University, Uxbridge, Middlesex, UB8 3PH, U K
e-mail: hamid.bahai@brunel.ac.uk
« 2 a parameter of the hepta-parameter spline distribution
«3 a parameter of the hepta-parameter sphne distribution
04 a parameter of the hepta-parameter spline distribution
day a number in the interval 0-365.25 w i t h step 0.125. 0 represents 0 am on January 1 ƒ„. transfer function
/ Y / ' ' ' ' the probability density function of random variable Xi
/?, the value of significant wave height observed at time
Hs significant wave height
Hs(t^ a random variable representing the Hs o f L performance (cost) function used i n network
training
nil a parameter of the hepta-parameter spline distribution
7772 a parameter of the hepta-parameter spline
distribution
77 number (size) of the data o artificial neuron output
input of the artificial neuron pdf probability density function R correlation coefficient R M S E root mean square error s artificial neuron bias
to the initial day for starting the T, simulation tl the day of the ith 3-hourly sea state
T- 3-hourly mean zero-up-crossing wave period T.(tD a random variable representing T, of
tZi a particular value o f the T, of the sea state o f t,-V output of the summing junction i n the
artificial neuron
J Mar Sci Teclmol (2007) 12:22-33 23
11' synaptic weight o f an artificial neuron Xj the output of the jth network, y = 1,. . . , 7 Xi TXQ given n(r,_,) = /z,_„ TXt,^,) = tz,_„
Hs(t,_,) = /;,_„ Hs{t,) = /?,, and = Z?,,, j'i ith simulated 71
Y mean of Zl ith observed r_ Z mean of z,
A a parameter of the hepta-parameter spline distribution
X a parameter of the hepta-parameter spline distribution.
1 Introduction
The knowledge of wind wave characterisdcs including the mean zero-up-crossing wave period (T^) is essential in a number of ocean engineering activities, such as those i n marine operations that call f o r real-dm'e or on-line forecasting o f the sea state characteristics. These characteristics have previously been simulated using analytical models by a number of researchers. A n alter-native simulation approach is based on artificial neural networks (ANNs), which have been applied effectively in predicting natural phenomena containing undeter-mined interrelationsftips among their physical parame-ters. Some recent apphcations of A N N s ' for predicting wave characteristics include the use of a three-layered feed-forward A N N to obtain the significant wave height (Hs) and 71 (output) f r o m wind speeds (input), resulting in an appropriately trained network which could pro-vide satisfactory results for certain types of predictions. Other types o f models for on-line prediction o f wave heights include^ first-order autoregressive moving aver-age ( A R M A ) and autoregressive integrated moving average ( A R I M A ) models. I t has generally been shown^ that A N N s result i n more accurate prediction o f wave heights than the two time-series models when shorter prediction intervals are involved; but for long-range predictions, both approaches appear to have similar performances. Other types of different neural network strategies^''' using various hidden layers have been em-ployed to forecast the sea state characteristics f o r various time intervals i n the future, and the suitabifity of A N N s has been demonstrated^ through verifying the short-term forecasts with observed data. However, the results of the simultaireous forecasts exhibited less ac-curacy than those obtained f r o m the separate predic-tions. I n nearly all applications of A N N s for predicting wave characteristics, some inputs such as wind speed or wave characteristics are given to various multilayered
networks to be trained for outputting wave character-istics such as or Hs.
A different approach has recently been adopted'' for the simulation and prediction of significant wave height using a set of feed-forward networks w i t h two hidden layers to estimate the parameters of the conditional dis-tribution of a desired future Hs given the immediate eight preceding Hs values to simulate the 3-hourly Hs sequences. The work reported in this article, unlike i n previous work, takes into account the stochastic nature of T- and introduces a conditional probability distribu-tion f o r 71 given some T. values and Hs values. A N N s are utilized to calculate the maximum likelihood esti-mates (MLEs) of the parameters o f t h e distribution f r o m the outputs of the networks trained to minimize the M L E function. Once the parameters are estimated, the simulation and prediction of are carried out through generating random variates f r o m the conditional distribution.
I n this study the A N N s have been trained with an optimization technique caUed simulated anneahng (SA), and the conditional distribution has been approximated by a proposed distribution called the hepta-parameter spline, described i n Appendix A . This distribution is utilized i n the prediction of the 71 o f the 3-houiiy sea states.
2 The data
The Hs and T, values used as input data i n this work were those of the 3-hourly sea states observed between 1978 and 1999 at US National Oceanographic Data Center Buoy 46005 in the Northeast Pacific (near 46°N 131°W), downloaded f r o m http://www.nodc.noaa.gov/ BUOY/46005.html. We expected to have 22 x 8 x 365.25 = 64284 observed Hs values and 64284 observed 71 val-ues; however, some are missing. The details of the miss-ing Hs and 71 values have been reported by Anderson et al.' The missing Hs or 71 values whose immediate pre-ceding and immediate following 3-hourly observed data were available, were interpolated linearly; the others re-main missing. I n total, 55922 values became available f o r training the networks. The training code has the option o f using all or a part of the observed data as the target f o r training. To increase the speed of training, nearly half the data, i.e., a vector of size 1 x 25 000, was used as the target vector during training the A N N s . This set of data (see Fig. 8) was observed f r o m Jamtary 1, 1978, to mid-1986. The rest o f t h e data, i.e., those ob-served f r o m inid-1986 to 1999 were used f o r testing the trained A N N s (see Fig. 9).
24 J Mar Sci Teclmol (2007) 12:22-33
3 Artificial neural networks
A N N s belong to that type o f information-processing system that allows an approximation to be made o f the nonlinear behavior (Haykin'*) that is the characteristic of geophysical processes (Makarynskyy et al.^). The basis of an A n n is the concept of neurons, as found i n bio-logical neural networks. Artificial neurons are extremely simple abstractions of biological neurons. A simplified mathematical model o f such a neuron is shown i n Fig. 1. The unit takes an argument v that is formed as the sum of a weighted input ( i i ' x and a bias; output o is produced by the means of a transfer (activation) func-tion denoted by which is typically a logsig funcfunc-tion, as described i n Appendix B. For the neuron illustrated in Fig. 1 we have:
0 = ƒ,,(!') = /;,.(wx/7 + ^ ) , (1)
where p is the input to the artificial neuron, w is the synaptic weight, s is the bias, ƒ , is the transfer function, V is the output of the summing junction, and o is the output of the artificial neuron.
A N N s are classified into several categories including multilayer feed-forward networks, which constitute one of the most common neural network modeling tech-niques. The current work uses this category of A N N . The details of A N N modeling are discussed in many references, e.g., Haykin.^
The weights and biases o f the networks are deter-mined through a process known as training. During this process the weights and biases are continually modified by means of a training algorithm to optimize a desired function called a performance (cost) function. I n this study, two training algorithms were tested: the gradient descent algorithm and simulated annealing. The stochas-tic algorithm of simulated annealing proved to be a much faster training algorithm for our purpose than the deterministic algorithm of gradient descent.
Fig. 1. A n artificial neuron with a single scalar input and a bias
• 0 Springer
4 Simulated annealing
Simulated annealing (SA) is a global optimization tech-nique for multivariable problems. Based on the princi-ples of thermodynamics, it is motivated by an analogy to annealing in solids. The idea stems f r o m the work of Metropolis et al.' Kirkpatrick et al.'° applied the idea to optiiuization problems and used SA to search for a fea-sible solution space and convergence to an optimal solution.
Figure 2 shows the flowchart of the standard SA al-gorithm used for training the networks in the present work. As in the process of physical anneahng, a param-eter called temperature is an important component of the simulated aimealing algorithm. The algorithm f o r a maximization problem starts at a high temperature f r o m a valid solution; new states are randomly generated for the problem, and the associated cost function, which is to be optimized, is calculated. A new state is randomly
Initialization (CurTent_solu[ion, temperature)
Calculation ofthe CurrentCost
LOOP
NewState
Calculation of tiie New_Cost
IF (Cuiient_Cost - New^Cost) < 0 THEN
Current_State= New_State
ELSE
Current Cosl - New Cost
IF Exp = =—>Random (0, 1)
THEN
— Accept
Current_ solution ~ New_ solution
ELSE
-- Reject
Decrease the temperature
EXIT when STOP CRITERION
END LOOP
J Mar Sci Teclmol (2007) 12:22-33 25
cliosen and "Current-Cost - New-Cost," i.e., tlie differ-ence in tfie cost function, is calculated. I f the differdiffer-ence is negative, the new state is accepted. This forces the system toward a state corresponding to a local, or possibly a global, maximum. However, most large optimization problems have many local maxima, and the optimization algorithm is hkely to become trapped in one of them. To get out of a local maximum, a decrease of the cost function is accepted based on a probability given in statistical thermodynamics, i.e.,
f Current _Cost - New _Cost"\ .
I I exp IS greater than a temperature
uniformly distributed random number i n the range 0 - 1 , the new state is accepted, despite its being a worse solu-tion. Using this approach, i t is most probable that all states of the system will be explored, including the global optimum. There are different methods for reducing the temperature iteratlvely in the cooling schedule. The method used i n this work, first suggested by Lundy and Mees," conducts only one iteration at each tem-perature and decreases the temtem-perature very slowly. I t should be pointed out that various researchers, including Burke and Kendall,'" have shown that the cost function can be responsible for a large proportion of the algo-rithm execution time. For an overview of SA, there are many references, including Dowsland'^ and Salamon etal.'" "
5 Specifications of the networks
The seven parameters of the hepta-parameter spline dis-tribution used for approximating the conditional distri-bution of 71 are estimated f r o m the outputs of seven trained networks, each providing one output for each set of inputs; the output of each network is used to estimate one parameter. Based on the experiences of Bazargan,' who successfully used feed-forward A N N s with two hid-den layers for simulating significant wave height, several different feed-forward network architectures, including
10 X 5 X 10 X 1, 10 X 5 X 5 X 1, 10 X 3 X 10 X 1, 10 X 1 X 10 X I , and 1 0 x 1 x 1 5 x 1 , v/ere tried and finally the size 10 X 1 X 15 X 1 (Fig. 3) was chosen for all o f the seven neural networks.
The T. values simulated using this architecture were closer to the observations. Each network, whose single
>|lW(l,ll[
r~bii)~i z
b[3)
output is used for estimating one of the parameters of the distribution, has ten inputs and two hidden layers having 1 and 15 neurons respectively. The output layer has one neuron. The transfer functions of the two hidden layers and the output layer were of the logsig, logsig, and purelin type, respectively. Networks having one hidden layer could have been used, but the hidden layer would have needed a larger number of neurons.
5.1 The target vector, input matrix, and output vector of the networks
As pointed out earlier, a vector of observed T, values having the size 1 x 25 000 was used as the target vector of each network. Each value in the vector is an observed mean zero-up-crossing period related to the sea state o f tj denoted by r.(f,), in which i.e. the day of the /th 3-houiiy 71 value satisfies Eq. 2:
;,. = ?o+0.125/ for / = 0,1,2,3,.., (2)
where t^ is a number in the interval 0-365.25 representing the starting day of the period whose 3-hourly sea states are to be predicted. The maximum value o f / i n this case is 25 000.
To increase the efficiency of training, the values of the target vector were scaled as follows:
rZsca,ed=0.8X r ^ - m i n ( 7 ; ;
m a x ( 7 : ) - m i n ( 7 ; ) -4-0.1 (3)
where min (71) and max ( T J are the minimum and the maximum o f the observed T- values, which are 3.7 s and 17.5 s, respectively.
The input data f o r training the networks consisted o f a matrix of size 10 x 25000. The matrix was created in such a way that f o r any 71 in the target vector, the f o l -lowing ten rows were i n the corresponding column o f t h e matrix:
Rows 1 and 2: r,(?,_2), T^tf-i), i.e., the tvv'O immediately preceding values, scaled by Eq. 3.
Row 3: Hs(ti_t), i.e., the immediately preceding Hs value.
Row 4: Hs{tD, i.e., the Hs o f the sea state of
Row 5: Hs(ti+t), i.e., the immediately following 3-hourly Hs value.
The significant wave heights were scaled as follows:
Fig. 3. The chosen artificial neural network ( A N N ) architecture
^^^sca,ed = 0.8X Hs-min^Hs)
max ( H s ) - min (77^) + 0.1 (4)
26 J Mar Sci Technol (2007) 12:22-33
where min(i75-) and max{Hs) are the minimum and the maximum of tlie observed Hs values, which are 0.2m and 13.6m, respectively.
The experiments of Bazargan'' showed that i f the time and season values of the desired T, are given to the net-works as a part of the input instead o f the day on which the value occurred, the convergence of the networks increases and their size decreases. Therefore, Rows 6 through 10 were added to the input:
Row 6: the time of the desired r_; this value lies i n the interval 0-24, but its scaled value, which lies i n the interval 0 - 1 , was used (see Appendix C) to increase the efficiency of training.
Rows 7 though 10: Four fuzzy values describing the season.
The season is usually considered as spring, summer, fall, and winter ( 1 , 2, 3, and 4); however, this work ex-presses it as a fuzzy value, such that a particular day o f the year belongs to all the four seasons with four mem-bership values, each lying i n the interval 0 - 1 . Therefore, four different membership values were assigned to each day of the year. For example the vector ( 1 , 0, 0, 0) was assigned to 7:30 am of the 127th day o f t h e year, or 5th May (day = 127.3125), and (0.0246, 0, 0, 0.9754) was assigned to 6 am on 7th February, or 6 am of the 38th day of the calendar (day = 38.25). The fuzzy diagram of the membership values is shown in Fig. 4. Appendix C shows how the membership values and the time are cal-culated f r o m the day.
The output o f each of the seven networks is a real number used to estimate the seven parameters o f the probability density function (pdf) of the spline distribu-tion, described i n Appendix A , for approximating the conditional p d f f o r / = 8 to 25000, where Xi repre-sents TXtd given
r,(?,_,) = te,_„ r,(?,,,) =
tz,._„ Hs{U_,) = /?,_„ Hs{t) = hi, and i75(?,.+,) = tó- is a particular values > a. a s u m m e r a u t u m n w i n t e r • day 4"' 5"'
Feb May Aug Nov Feb May
Fig. 4. A fuzzy diagram for season membership values
of T^t^, Hs(tD is a random variable representing the Hs of and hi is a particular value of HsQ^. EachfJ"^ be-longs to a T ; in the target vector.
Since the training algorithm (in our case SA) requires a cost function to be defined f o r optimization, the func-tion, based on the function of the maximum likelihood estimation ( M L E ) method, given by Eq. (5) was selected f o r this purpose:
L= t l°gA'/'""''- (5)
6 Training process
Because L i s a function of the p d f s parameters estimated f r o m the outputs of the networks, and the outputs of the networks depend on the networks' weights and biases, our objective during the process of training was to change the weights and biases continually i n order to maximize L. However, to get more accurate resuhs, the root mean square error (RMSE) of the observed and predicted values f r o m the corresponding distributions was simultaneously monitored using Eq. 6 f o r a period starting f r o m January 1, 1978.
The weights and biases o f the seven networks were changed randomly during the training process. For each 71(0 in the target vector, each of the seven networks would create an output, i.e., seven outputs for each T^tl). The seven values were used i n estimating the seven parameters o f f r o m Eqs. 13-20 of Appendix B. The new value o f the cost function and that of the R M S E were then calculated. The new weights and biases were accepted or rejected according to the following rules: — I f the new values of the cost function and the R M S E
have improved in comparison to the previous itera-tion, the new weights and biases are accepted. — I f the cost function and the R M S E have both become
worse, the increase in R M S E is considered. I f it has not exceeded a specified small value, the new weights and biases are accepted w i t h the probability men-tioned in the SA algorithm, but i f the R M S E has exceeded the specified small value, the weights and biases are rejected.
— I f the cost function has improved but the R M S E has become worse by a small amount, the weights and biases are accepted. I f the cost function has improved but the R M S E has become considerably worse, the weights and biases are rejected.
— I f the R M S E has improved but the cost function has become worse, we accept the weights and biases with some probability. I n other words we first check the
J Mar Sci Technol (2007) 12:22-33 27
R M S E value; i f it has become considerably worse, the new weights and biases are rejected no matter whether the cost has improved or not. But i f the R M S E has become worse by a small value or has improved, then we look at the cost and accept the biases and weights i f the cost has improved. They are accepted w i t h some probability i f the cost has be-come worse.
I n order to prevent over-fitting, the simulated T_ val-ues f o r the period 1978-1999 were occasionally plotted versus the corresponding observed values and the plot was checked visually. Over-fitting is a situation that may occur during network training in which the error on the training set was driven to a very small value, but when some new data is presented to the network, the error becomes considerable. The process of changing weights and biases continues until certain measures f o r terminat-ing the trainterminat-ing are satisfied. Our measure of termination was to train the networks in such a luanner that L and the R M S E discussed earlier were optimized. When train-ing was terminated, the cost function had reached 70529 and the R M S E f o r the period whose data were used dur-ing the traindur-ing was 0.727 s.
7 Simulation of T.
To simulate the T- values of the sea states of a future period, as wefi as the following data need to be given as input to the M A T L A B code:
1. TJyt^.^), i.e., the 71 o f t h e sea state 3h before 2. T-iio-i), i-e-, the 71 of the sea state 6h before /„, 3. Hs(tg_^), i.e., the Hs of the sea state 3 h before /„, 4. and 5. Hs(ta), Hs{to+i), i.e., the Hs values o f the sea
states o f t-Q and 3h after tg, respectively.
Hsito) and Hsitg^^) will take place in the future, but they are forecast using the conditional distributions of Hs{t^ and Hs(t(,^i), whose parameters are readily given by feed-foreword A N N s o f size 13 x 24 x 13 x 1, already trained f o r simulating Hs by Bazargan.''
The f o u r fuzzy membership values and the time are calculated f r o m using the instructions given in Appendix C. The outputs of the trained networks are obtained and the parameters of the conditional distri-bution of the desired wiU then be estimated by a M A T L A B code which uses the equations given in Appendix B. The conditional distribution is now completely known.
Throughout this work, it is assumed that the condi-tional distributions of 77^ and do not change f r o m one year to another; i.e., the conditional distribution o f is the same as that of TXti + 365.25).
The inversion method (Naylor et al."') was used to simulate the values o f the sea states of the desired pe-riod one after another f r o m the hepta-parameter spline distributions approximating the conditional distribution of the values. Figure 5 shows the 22-year observed T, values (gray) and the results of a typical simulation o f T- values (black) f o r the Northeast Pacific starting f r o m /„ = 00 am January 1 1978, to 1999, given the necessary data. Figures 6 and 7 show the observed and the simu-lated values separately.
The traiiung set o f data and the test data, together with their corresponding simulated data, are shown in Figs. 8 and 9. As the figures show, the model has learned to imitate the behavior of the observed values.
Figure 10, showing the cumulative distribution func-tion of the above observed and the simulated values, conveys the degree o f agreement between the observed and simulated T, values. I t can be shown easily that
_10 Ll 8
* "t/'^'^ '"'''r.r-i ^ ' y W
0 1000 2000 2000 4000 flOOü 6000 7000 8000 Titrietpayj
Fig. 5. The observed (gray) and simulated (black) r . s of a region in the Pacific from 1978 to 1999 using ANNs. The gray data points have covered the black ones
0 1000 2000 3000 4000 5000 6000 7000 6000 Timet Dav)
Fig. 6. The T, values observed from 1978 to 1999
28 J Mar Sci Teclmol (2007) 12:22-33
A typical simulation for 1976-99
: : QOO^A^^r^JfjOO 6000 7000 8000
Fig. 7. The T. values simulated from 1978 to 1999
16
T r a i n i n g d a t a li^jid
14 f- the c«rrespondini:j simulated Tz's
. ^ I
f :
- I I I
0 500 1000 1500 2000 2500 3000 3500 4000
Time (dav)
Fig. 8. Training data observed from 1978 to 1986 (black) and the corresponding simulated data (gray)
Fig. 10. The empirical cumulative distribution function of the observed and simulated T. values from 1978 to 1999
Iraving simulated the 71 value of a future sea state a suf-ficiently large number of times, an appropriate predic-tion f o r T. is the mean of the simulated values. This is illustrated i n Sect. 9.
8 Network efficiency
The efficiency of the artificial neural networks was cal-culated in terms of three indices, i.e., the RIVISE, the correlation coefficient (7?), and the scatter index (SI), according to Eqs. 6-8, respectively:
R M S E : (6) 16 14 12 -•10 3 6 4
Test data and
the corressponding :5irnulated T j ' s
1^ 15
3.500 4000 4500 5000 SSOÜ 6000 8500 7000 7500 8000
Tirne(day)
Fig. 9. Test data observed from mid-1986 to 1999 (black) and the corresponding simulated data (gray)
R- Sl-J^(z,^Z){y-Y) i=\ X ( z , - Z ) ^ t ( T , - F ) ^ 1=1 1=1 R M S E (7) (8)
where z, is the /th observed T., j ' , is the /th simulated T., n is the number (size) of data, Z is the mean of z„ and
Y is the mean of j ' , .
For the 22-year typical simulation shown i n Fig. 5, these indices were: R S M E = 1.4 s, 7? = 0.5, and SI = 0.21.
A n illustration of short-term forecasting is now given. Since the parameters of the conditional hepta-parameter spline distribution o f each of the T. values in the target vector (which are the observed T, values of all of the
J Mar Sci Teclmol (2007) 12:22-33 29
3-hourly sea states f r o m 1978 to mid-1986) could readily be calculated f r o m the trained outputs of the networks, the above T, values were forecast i n the following three ways f r o m the corresponding conditional hepta-parameter spline distribution, given the necessary data: 1. The mode (most probable value) of its corresponding
spline distribution
2. the mean of the distribution, and
3. a variate generated directly f r o m the distribution. Table 1 shows the results of the calculations of the three indices o f efficiency f o r these methods using the predict-ed values and the observpredict-ed values. According to the indices in Table 1, the mean of the distribution is a better forecast.
The scatter plot of the observed values and the predicted values obtained f r o m the means of the distri-butions is shown i n Fig. 11.
The quantile-quantile (QQ) plot in Fig. 12 shows that the samples of the observed and predicted T, values could be considered f r o m the same distribution type.
As another typical illustration. Fig. 13 shows the ob-served and simulated data f o r 1978. Since the obob-served
Table 1. The efficiency indices for three different forecasts of the immediately subsequent T. for 1978-1986
Forecast RMSE (s) R Sl
Most probable 0.87 0.86 0.12
Expected value of T. 0.72 0.88 0.10 The two variates generated 1.12 0.67 0.16
1.11 0.69 0.16
RMSE, route mean square error; R, correlation coefficient; SI, scatter index
values have not been reported f o r a period of the year, no T, values were generated for that period, and a gap is seen in the figure.
9 Prediction of sea states
The above procedure, together with the work of Bazargan,'' can simulate and predict the characteristics of the sea states at 3h, 6h, 1 day, or 1 week, etc. f r o m an initial dme. As an illustration, suppose the sea state
14
12
10
t
-UbsGtved arid Predicted Tz sf(:ir197B-mid 1986, QQ Plot
8 10
Obsetved T i (s)
12 14
Fig. 12. The 3-hourly T, values observed from 1978 to mid-1986 versus the predicted 71 values, quantile-quantile plot
Obsetved and Predicted 3-hourly T i for197S-mid 1936
6 s 10
Ob.served T i (s)
14 16
Fig. 11. The scatter plot o f t h e observed 3-hourly T. values in the period 1978 through mid-1986 versus the predicted T- values ob-tained from the mean of the distributions
14
One-year isirnulation of T i
100 200
TimeCDay)
300
Fig. 13. The Tz values for 1978; the observed Tz values are in gray and the simulated values are in black
30 J Mar Sci Teclmol (2007) 12:22-33
Table 2. Inputs for the sea state prediction
Time Date Hs{m) r.(s) 12 pm Jan. 1 3.1 6.7 9 pm Jan. 1 3.2 6.8 6 pm Jan. 1 3.0 3 pm Jan. 1 2.6 12 n Jan. 1 2.7 9 am Jan. 1 2.0 6 am Jan. 1 1.7 3 am Jan. 1 1.8
Hs, significant wave height; T., 3-hourly mean zero-up-crossing wave period
Table 4. The results of generating the sea state ten times for 0 am on December 31 Run no. Hs(m.) 1 4.2 8.4 2 9.3 10.4 3 4.7 7.9 4 4.3 7.6 5 2.7 6.8 6 4.2 7.5 7 1.1 6.7 8 6.9 8.0 9 9.1 9.1 10 5.8 10
Table 3. The results of generating the sea state ten times for 12
noon on January 2 Run No. Hs(m) r,(s) 1 2.5 7.6 2 3.1 7.9 3 3.7 7.4 4 4.3 7.9 5 4.1 7.0 6 3.8 8.3 7 5.4 7.3 8 2.8 6.6 9 3.9 8.3 10 5.5 7.4
of day = 1.5, i.e., 12 noon on January 2nd, for the region in Northeast Pacific is to be predicted. Table 2 shows the observed Hs and T. values, all to be given as input to the M A T L A B code used f o r the prediction. A f t e r inputting the data into the code that calculates the parameters of the spline distributions f r o m the outputs of the trained networks and generates the Hs and T, values as random variates f r o m the corresponding distributions, the sea states on January 2nd f r o m 3 am through 12 nooit were generated ten times by running the code. The results of the simulation for 12 noon on Jaiutary 2nd are shown in Table 3. For 12 noon on January 2nd, the luean of the ten sets of generated data for Hs and T, are 3.91 m and 7.57 s respectively. These results compare well with the corresponding observed values, because the mean of the 22 values o{ Hs and T, observed at 12 noon on January 2nd for 1978-1999 are 3.89m and 7.67s, respectively. For practical purposes, the number of runs would be increased to gain more accuracy.
As a second illustradon, after giving the same input as i n the former illustration, the sea state for 0 am on December 31 was simulated ten times as shown i n Table 4. A suitable prediction f o r the sea state given the Hs values of this time of the year could be the sea state characterized by the mean of the Hs values and the T,_
values of Table 4, respectively, i.e., Hs = 5.2m and -8.2 s. The mean of the observed values of Hs and for 0 am on December 31 for the period 1978-1999 were 3.46m and 7.52 s, respectively. Comparing the first and the latter illustration, it is concluded that the closer i n time the input data and the predicted value are, the more accurate the prediction.
10 Joint distribution of Hs and
A consistent method of choosing the height and period for the adopted design wave may be crucial for reliable designs of offshore structures. W i t h this objective, the joint pdf of the height and period is considered (Haver"). As an illustration, suppose the joint distribution o f Hs and T. of the sea state at 6 am on A p r i l 1, 2004, given the Hs values and 71 values of January 1, 1978, is to be developed. A two-dimensional histogram of relative fre-quency of a sample of Hs(t^ and T-(ti) approaches their joint distribution i f the sample size is large enough. The Hs values and values of the period 0 am on Jaiuiary 2, 1978, through 6 am on A p r i l 1, 2004, could be simu-lated using the networks trained for Hs and 71. The simulated data is a sample f r o m the conditional joint distribution of:
Hs (0 am Jan2,\978), Hs (3 am Jan 2,1978),..., Hs {6 am A p r i l 1,2004),
" T (0 am Jan 2,1978), T, (3 am Jan 2,1978),..., X (6 am A p r i l 1,2004),"
given the 3-hourly observed Hs values and r„ values for January 1, 1978.
Indeed the simulated sea state characteristics (Hs and T^) f o r 6 am on A p r i l 1, 2004, are random variates f r o m the conditional joint distribution o f Hs and T, related to 6 am on A p r i l 1, 2004, given the observed sea states for
J Mar Sci Teclmol (2007) 12:22-33 31
January 1, 1978. However, the difference between these two times is rather long; therefore, i t is quite logical to neglect the information conveyed by the sea states o f the starting day and assume that the variates come f r o m the joint distribution of the characteristics o f the sea state at 6 am on January 1, 2004, rather than coming f r o m the conditional joint distribution o f the characteristics.
By repeating the above simulation a sufficiently large mimber o f times, a largish sample o f the characteristics of the sea state at 6 am on A p r i l 1,2004, can be simulated and f r o i u which a two-dimensional histogram o f relative frequency that represents approximately their j o i n t dis-tribution can be developed.
11 Conclusions and suggestions
The above two properties imply that there are two points, say a and e {a < e) for which the value o f f ( x ) approaches zero f o r A' < a and x > e and also there is a point c i n the interval <a-e> f o r w h i c h / ' ( c ) = 0 (e.g., the maximum of the function). There is no a priori information about the shape o f the distribution; therefore, a parametric density function must be found that approximates a number of probability distributions as its parameters are changed. The cubic interpolations on the four subintervals of Fig. 14 were assumed to be a good approximation for f(x). Note that a, b, c, d, e, /»,, and nh were chosen such
that:
b = ^ , d = ^ , m , ^ f ' { b ) , and m,=f'{d).
1. Multilayer feed-forward A N N s could be used effi-ciently in estimating the parameters o f a conditional distribution used to simulate the mean zero-up-crossing wave period o f the sea state of a desired time.
2. The proposed hepta-parameter spline distribution approximates well the conditional distribution used for simulating the T. values o f Northeast Pacific 3-hourly sea states.
3. The mean o f the simulated values obtained f o r 71 values in several runs is a good forecast f o r the 71 o f sea states not too far i n the future.
It should be noted that networks o f size 10 x 1 x 15 x 1 are not necessarily the best network architecture that could be trained f o r this purpose. More experiments with other network structures, other inputs, other mea-sures o f terminating the training process, and other probability distributions, along with performing sensi-tivity analysis, may result in better-trained networks. Other algorithms such as genetic algorithms could also be tried instead o f A N N s .
Acknowledgments. The helps of Mr. Amin Aminzadeh and M r . Ali Bazargan is highly appreciated.
Appendix A
Hepta-parameter spline distribution
For the simulation o f
I'-C;).
the distribution o f TJj^ given r,(/,_,)?z,-„ Wi_^ = tz,._„ Hs(ti_,) = A,^,, 775(0 = Hs{ti+t) = /7,-+,, whose pdf is denoted hyj{x) was used.y(x) should satisfy the following two properties:1. l i m / ( x ) = 0
2. Smoothness.
Figure 14 also shows «i, a^, « 3 , a^, A, A, aud 7772, the
parameters o f the function used f o r the approximation of the conditional pdf.
The function g(x), the proposed p d f f o r approximat-ing/(A-), is defined as follows, considering Fig. 14:
0 Cubic-Int and Cubic-Int and Cubic-Int. and Cubic-In and with slope = 0 .oJ with slope 777, > 0 with slope 771, > 0 with slope = 0 with slope = 0 with slope 7772 < 0 with slope 77 72 < 0 x<a a<x<b '<X<C •<x<d with slope = 0 0 d <x<e x>e where a, b, c, d, e, c^, ci^, Ö 4 , 7 7 7 , , are illustrated in Fig.
14, and Cubic-Int. stands f o r the cubic interpolation of each pair o f the points mentioned above according to the spline method.
Having two points A-, and A ' , and the values o f a func-tion/(A-) and its derivative at A , and A2, the formula f o r approximating/^A'^ according to the spline method is given by:
32
J Mar Sci Technol (2007) 12:22-33
f(b) = 32
./
slope = m. )>a3 = f(c) f(d) = a4 a b e ai ai+XA/2 ai+X<
A,A
•U—
ai+(0.5+0.5?i) A (1-^)A — e ai+AFig. 14. The intervals on which the conditional distribution of T IS defined and the parameters of the distribution
rW = / ( A - 2 ) X - X , ( x - x , f ( x - X 2 ) 3 / ( x 2 X , - X , x - x , X , - X , ( x , - x , )
+ f / ' ( , - J _ M f L ) l ( f Z f 2 ) l ( £ z f J
^ ' ^^-^-2) ( x , - x , f Since J .^C-^')"^-^' = 1,, it can be shown that:X^a^ ( l - A ) A a
(9)
2 ' ^ - t ^ - • = 1 (10)
Hence, one of the parameters a,, ci„ a„ X, and A can be calculated f r o m the other four, and the mnnber o f inde-pendent parameters is reduced to seven; thus ^(x) has seven independent parameters and that is why the name hepta-parameter spline was chosen f o r it. Arbitrarily, aj was calculated f r o m a,, a„ X, and A. Finally, since all the interpolated values of ^(x) have to be positive, it can be shown that the following two inequalities must also be satisfied: 0 < m ,
A
' AA and -6a. (11) (12)The distribution has the capability of approximating many o f the conventional distributions, including the chi-squared, the positive normal, and the three-parameter Weibull distributions. However, an attempt to use the three-parameter Weibull distribution instead of the hepta-parameter distribution did not yield good results.
Appendix B
Calculation of the distribution parameters f r o m the networks outputs
Let Xj denote the output of the jth network for 7 = 1, 2, . . . , 7. fl| is positive, therefore we cannot use the outputs of any of the seven networks f o r calculating a, directly because the range of the outputs o f t h e networks is (-00 to °o) because they have the transfer function purelin f o r their output layer. Hence a transformation was applied to x, to arrive at a positive immber: = x , l The relationships used to calculate the parameters f r o m the outputs o f the networks are given i n Eas 13-20. Constraint a, > 0 0 < (72 < « 3 Relationship « 1 = a, = Ö3 logsig X2 A > 0 0 < r/4 < (73 0 < „ 7 , A ' AA A > X32 « 4 = 03 logsig X4 6(72 , ' " 1 = ^ l o g s i g X5 - 6 ( 3 4 < 7 ? 7 , < 0 777, = - -60,, - logsig Xg ( l - ^ ) A ^ ( l - A ) A 0 < A < 1 A = logsig X7
F r o m Eqs. 9, 13, 14, 15, and 18 we have:
(13) (14) (15) (16) (17) (18) (19)
2x3-(logsigx, )(logsigx2) + 0.5 + (1 - logsigx^ )(logsigx4) (20) where
logsigx^.
l+e - , j = l,...,l. (21)
I n practice, i n the first step is calculated by Eq. 20, then Ö2 and ,r/4 are calculated f r o m Eqs. 14 and 16. Other parameters are calculated simply using the appropriate relationship.
J Mar Sci Technol (2007) 12:22-33 33
Appendix C
Calculation of the time and four season values f r o m day values
The time at which the sea state being predicted occurs is calculated f r o m the following M A T L A B command: time = day-floor (day)
where day is a mimber in the interval 0-365.25, with step 0. 125. representing the day on which 3-hourly sea states occur; day = 0 represents 0 am on January 1.
The f o u r fuzzy membership values related to the sea-son i n which the sea state occurs are calculated using the foUowing four M A T L A B instructions:
1. springvalue = Trianglefun(abs(day - midspring)/ (365.25/4))
2. sunmiervalue - Trianglefun(abs(day - midsummer)/ (365.25/4))
3. fallvalue = Trianglelfun(abs(day - midfall)/ (365.25/4))
4. wintervalue = Trianglefun(abs(day - midwinter)/ (365.25/4))
where
Trianglefun is the following M A T L A B function: f u n c t i o n v = Trianglefun(day);
V = (1 - abs(day)) (abs(day) < 1)
+ (1 - abs(day - 4)) (abs(day - 4) < 1)
and, since the 36th day of the year (5th February) is midwinter,
midwinter = 36; and also:
midspring = 36 + (365.25/4) = 127.3125 (5th May) midsummer = 36 + (2 x 365.25/4) = 218.6250 (6th August) m i d f a l l = 36 + (3 X 365.25/4) = 309.9375 (6th November) References
1. Deo M C , Jha A , Chaphekar A, Ravikant IC (2001) Neural networks for wave forecasting. Ocean Eng 28:889-898 2. Agrawal JD, Deo M C (2002) On-line wave prediction. Mar
Struct 15:57-74
3. Makarynskyy O, Pires-Silva A A , Makarynska D , et al (2002) Artificial neural networks in the forecasting of wave parameters. 7th International Workshop on Wave Hind-casting and ForeHind-casting. Banff, Alberta, Canada, Oct 2002, pp 514—522. http://www.waveworkshop.org/7thwaves/papers/ makarynskyy_etal.pdf
4. Makarynskyy O (2004) Improving wave prediction with arti-ficial neural nets. Ocean Eng 31(5-6):709-724
5. Makarynskyy O (2005) Artificial neural networks for wave tracking, retrieval and prediction. Pacific Oceanogr 3(1): 21-30
6. Bazargan H (2007) Neural network-based simulation of sig-nificant wave height. Part of a PhD thesis (http://hdl.handle. net/2438/379). School of Engineering and Design, Brunei University, U K . In: Proceedings of 26th international confer-ence on offshore mechanics and arctic engineering, June 10¬ 15, 2007, San Diego, USA
7. Anderson CW, Carter D , Cotton PD (2001) Wave climate variability and impact on offshore design extremes. Report for Shell International and the Organisation of Oil and Gas Producers. http://info.ogp.org.uk/metocean/JIPweek/ WCEReport_2sided.pdf
8. Haykin S (1999) Neural networks: a comprehensive founda-tion. Prentice-Hall, Upper Saddle River, NJ
9. Metropolis N , Rosenbluth A W , Rosenbluth M N , et al (1953) Equation of state calculation by fast computing machines. J Chem Phys 21:1087-1092
10. Kirkpatrick S, Gelatt CD, Vecchi MP (1983) Optunization by simulated annealing. Science 220(4598):671-680
11. Lundy M , Mees A (1986) Convergence of an annealing algo-rithm. Math Prog 34:111-124
12. Burke E K , Kendall G (1999) Evaluation of two-dimensional bin packing problem using the no fit polygon. In: Proceedings of the 26th International Conference on Computers and Industrial Engineering. Melbourne, Australia, 15-17 Dec 1999, pp 286-291
13. Dowsland, K A (1995) Chapter 2. Simulated annealing. I n : Reeves CR (ed) Modern heuristic techniques for combinato-rial problems. McGraw-Hill, Berkshire
14. Salamon P, Sibani P, Frost R (2002) Facts, conjectures and improvements for simulated annealing. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA 15. Moins S (2002) Implementation of a simulated annealing
al-gorithm for M A T L A B . A Report to the Electronic Systems Department of Linköping Institute of Technology, Linköping Sweden, http://www.ep.liu.se/exjobb/isy/2002/3339/
16. Naylor H N , Balintfy JL, Burdick DS, et al (1968) Computer simulation techniques. Wiley, New York
17. Haver S (1987) On the joint distribution of heights and periods of sea waves. Ocean Eng 14(5):359-376
