Artificial Neural Network Model for the Evaluation of Added Resistance of Container Ships in Head Waves

Abstract

The decrease in ship added resistance in waves fits into both the technical and operational measures proposed by the IMO to reduce the emissions of harmful gases from ships. Namely, the added resistance in waves causes an increase in fuel consumption and the emission of harmful gases in order for the ship to maintain the design speed, especially in more severe sea states. For this reason, it is very important to estimate the added resistance in waves with sufficient accuracy in the preliminary design phase. In this paper, the possibility of applying an ANN to evaluate added resistance in waves at the different sea states that the ship will encounter during navigation is investigated. A numerical model, based on the results of hydrodynamic calculations in head waves, and ANN is developed. The model can estimate the added resistance of container ships with sufficient accuracy, based on the ship characteristics, sailing speed, and the sea state using two wave energy spectra.

1. Introduction

The evaluation of the ship added resistance in waves has increased in importance, especially from an economic and environmental protection point of view. Namely, added resistance in waves causes either a speed reduction or an increase in the required power and, consequently the fuel consumption, of a ship, which has a negative impact on its CO₂ emissions. The International Maritime Organization (IMO) has subjected the emission of harmful gases to increasingly stringent regulations with the aim of increasing the energy efficiency of ships and reducing CO₂ emissions. Considering the regulations adopted by IMO [1], the decrease of ship added resistance in waves could fit into both the technical measures for newly built ships and operational measures for existing ships in service. For this reason, it is very important to predict the increase in ship resistance due to waves to ensure the sufficient power required for sailing in actual sea states, but also to estimate the fuel consumption with the aim of reducing operating costs and emissions of harmful gases. Even though a ship is typically designed for calm water conditions, added resistance in waves should be considered as well in the ship design phase. When the ship is sailing in real conditions, the increase in required power to maintain the design speed, corresponding to calm waters, is traditionally accounted for through sea margin, which accounts for 15% to 25% [2] of the power required in calm water. However, that is a rough estimation based on the experience. For example, the sea margin due to wind and waves for an S175 container ship at the representative sea conditions, characterized by Beaufort number 6 for different wind speeds, was evaluated in [3], and it ranged from about 17% to 35%. The authors concluded that, for lower operating speeds, the increase in speed reduction, mainly due to head waves, was more pronounced. In addition, the evaluated sea margin was much larger when the initial operating speed was lower, accounting for the required increase in power up to 35%, despite the fact that the absolute value of the required power was obviously lower for lower speeds. Youngjun et al. [4] evaluated the actual sea margin of an LNG carrier for a given sailing route under wind and waves, by taking the fouling effects into account. The sea margin, based on the re-docking period, ranged from approximately 13% to 32% for the fouling effects that occurred within the period from one to five years. Degiuli et al. [5] demonstrated that the increase in the total resistance in waves, depending on the sea state, could be up to 23% for the design speed of 24 knots, and up to 36% for the slow steaming speed for the sailing route passing through the North Atlantic. These increases would be lower in the Mediterranean Sea due to lower sea states. Thus, Degiuli et al. [6] showed that the increase in the total resistance in waves for different sea states can be up to 5.5% for the design speed of 22 knots and up to 4.5% for the slow steaming speed.

Whether they are used to determine the power required for ship to maintain speed in real sailing conditions or the involuntary speed reduction, numerical calculations or model tests should be conducted. On the other hand, analytical expressions are more robust, but often of limited accuracy. Added resistance is considered a non-viscous phenomenon, which allows the application of methods and solvers based on the potential flow theory, but also allows an extrapolation of the results from model scale to ship scale without the scale effects. Indeed, the viscous part of added resistance can take a notable portion in added resistance only in short waves, when it can account for up to 20% of the added resistance for the model scale [7], and the viscous effects are even less present in full-scale. Namely, in short waves the portion of frictional resistance in total resistance in waves is larger than that of pressure [8]. Based on extensive study, Bunnik et al. [9] concluded that as long as strongly non-linear effects are not present, computational fluid dynamics (CFD), based on the viscous flow theory, does not contribute significantly to the accuracy of the numerical results of ship motions and added resistance in waves. On the other hand, CFD based on the viscous flow theory in comparison to the methods based on the potential flow theory increased the required computational time considerably.

As the IMO emphasized the need for the development of the so-called level 1 methods for the estimation of the required power in severe weather conditions [10], a relatively simple method for the determination of added resistance in waves was developed in [11,12]. The proposed empirical method requires the wave characteristics as the input, along with the main particulars determined based on the assessment of the influence of the ship’s main particulars and mass characteristics on added resistance in waves. The authors concluded that the peak position of added resistance was closely related to the heave natural frequency in long waves and that the amplitude of added resistance highly depended on the gyration radius. This finding was in accordance with [13], where the authors investigated the effect of ship characteristics on added resistance in regular waves and actual sea states. The determination of the added resistance in waves, as the time averaged second order wave force, required rather complex hydrodynamic calculations to ensure an acceptable accuracy of the results. Therefore, a model that allows for a simple, but sufficiently accurate and reliable evaluation of the ship’s added resistance while sailing at an actual sea state could be of great benefit, especially during the ship design phase. Bearing in mind that artificial neural networks (ANN) are successfully applied in various fields due to their exceptional ability to predict and classify, such a model could benefit from the advantages of ANN. The advantage of ANN over some other statistical tools lay in their ability to generalize, but also their ability to work with incomplete and big data. As they are very adaptable and have the ability to learn from examples, as in many technical fields, ANN are used in ship hydrodynamics to solve various problems with great success.

Gkerekos et al. [14] combined multivariable regression analysis and ANN, and established a model to predict ship fuel consumption for draught, speed, route length, and environmental conditions as input variables. The authors expected that the proposed methodology would enable the control of harmful gas emissions via the optimization of the ship operability. Based on the collected data from a VLCC (very large crude carrier), and using a multivariable regression analysis and ANN, a model for predicting the required ship power and fuel consumption in real time was developed in [15]. Thus, as part of the decision support system, this developed model could play a key role in increasing ship energy efficiency during navigation. A decision support system based on ANN and using sea trial data was proposed in [16]. The developed model is able to predict fuel consumption and route duration using the shaft speed and propeller pitch as input parameters.

Kim et al. [17] proposed an ANN model for the evaluation of ice resistance of various ice-going ships for different ice thickness. Im and Nguyen [18] proposed a control and management system based on a neural network for ship berthing. The main benefit of the proposed model is that there is no need to re-train the neural network for unknown ports that are topologically similar to the original port for which the training process was carried out. A similar model for docking an autonomous ship was developed in [19] to determine the required rudder deflection angle and propeller speed during docking. The authors showed that the developed model provided reliable results for constant and zero wind loads, while further improvement of the model was required in the case of dynamic wind loads. Prpić-Oršić et al. [20] developed a methodology to estimate the wind loads acting on the superstructure of container ships based on the results of a CFD analysis and general regression neural network (GRNN). Awad et al. [21] integrated an advanced neural network model with error back propagation in the so-called direct inverse neural network control (DINNC) system, in which the neural network can produce a signal to the dynamic system to cancel or dampen the roll motions caused by external disturbances. Cepowski [22] analysed different topologies of static feedforward neural networks and applied various learning algorithms to determine the coefficient of ship added resistance in regular waves. The neural network learned based on the experimental data of 14 models of different ship types. The author showed that the model, which was based on the multi-layer perceptron (MLP) network, could be of great use in the preliminary design of ships.

In this paper, a numerical model based on ANN was developed for the evaluation of ship added resistance in head waves. The proposed numerical model could have a practical benefit during the ship design phase. The results of the hydrodynamic calculations of added resistance in waves for various hull forms of container ship at different sea states and ANN form the basis of the developed model. Namely, ANN has the ability to learn from examples and to estimate relationships between the input data and solutions to the nonlinear multivariable regression problems. Hydrodynamic calculations are conducted using the boundary integral equations method (BIEM) based on the potential flow theory. The obtained numerical results were validated against the available experimental data and verified (i.e., the numerical uncertainty is assessed). Finally, the applicability of ANN in solving complex nonlinear problems in ship hydrodynamics is demonstrated.

2. Numerical Calculation of Added Resistance

Ship added resistance in waves is determined based on the potential flow theory under the assumption of incompressible, inviscid, and irrotational flow with fluid motions represented by the velocity potential $\mathrm{\Phi }$, satisfying the Laplace equation in the entire computational domain. According to Green’s second identity, the velocity potential at any point of the computational domain can be determined based on the distribution of singularities over the domain boundaries, and the gradient of the velocity potential corresponds to the fluid velocity. Within the applied 3D panel method, the ship hull is discretized using quadrilateral panels with the constant distribution of sources on each panel. This results in integral equations, which are then iteratively solved to determine the unknown strength of the distributed sources [23], velocity potential, and pressure on each panel by the modified Bernoulli equation. The normal component of velocity potential at the collocation points corresponds to the normal velocity. As the fundamental solution of the Laplace equation, the Kelvin type Green function that satisfies the linearized boundary condition on the free surface, the boundary condition on the seabed, and the radiation condition is utilized. Since the Green function satisfies the boundary condition on the free surface in the entire computational domain below the linearized free surface, the complementary interior hull problem is implicitly solved at the same time as the exterior one. For that reason, at the so-called irregular frequencies, an unphysical solution can occur [24]. In this paper, to shift the irregular frequencies towards higher frequency range, an extended BIEM, that imposes Neumann’s “rigid lid” condition, was applied. The damping term is introduced in the momentum equation through fictitious force, dependent on fluid velocity, in order to artificially include the energy dissipation due to fluid viscosity. By developing the velocity potential above the linearized free surface, and neglecting the square terms, the boundary value problem (BVP) of the first order is obtained:

$${\nabla }^2\mathrm{\Phi }=0,$$

(1)

$$g{\mathrm{\Phi }}_z+{\mathrm{\Phi }}_{tt}+\mu {\mathrm{\Phi }}_t=0\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{for} z=0,$$

(2)

$${\mathrm{\Phi }}_n=X_t·n\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{on} \mathrm{wetted} \mathrm{surface},$$

(3)

$${\mathrm{\Phi }}_z=0\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{for} z=-h,$$

(4)

where $X_t$ is the velocity vector, $n$ is the normal vector, index $z$ denotes space derivative, indices $t$ and $tt$ time derivatives, $\mu$ is small positive parameter that introduces the energy dissipation into the momentum equation, and ${\mathrm{\Phi }}_0$ is the velocity potential of the incoming wave. Additionally, velocity potential satisfies the radiation condition on a cylindrical surface at infinity. It should be noted that the Green function satisfies all boundary conditions, except the one on the wetted surface. In order to account for the forward speed, the Green function associated with the encounter frequency was used.

The velocity potential and Green function can be expressed as $\mathrm{\Phi }\left(P,t\right)=\mathrm{Re}\left\{\varphi \left(P\right){\mathrm{e}}^{-\mathrm{i}\omega t}\right\}$ and $\mathcal{G}\left(P,Q,t\right)=\mathrm{Re}\left\{G\left(P,Q\right){\mathrm{e}}^{-\mathrm{i}\omega t}\right\}$, respectively. Thus, the Green function represents the field of velocity potential at the point $P$ created by the source of unit density located at $Q$. The velocity potential consists of the incoming wave potential, radiation potential ${\mathrm{\Phi }}_j$, and diffraction potential ${\mathrm{\Phi }}_7$, as follows:

$$\varphi =-\mathrm{i}\omega {\displaystyle \sum _{j=1}^6{\zeta }_{aj}}{\varphi }_j+{\zeta }_a\left({\varphi }_0+{\varphi }_7\right),$$

(5)

where ${\zeta }_{aj}$ and ${\zeta }_a$ denote amplitudes of the ship motions and incoming wave, respectively. The source distribution, $\sigma ={\varphi }_n-{\varphi }_n^{\prime }$, can be determined by satisfying the boundary condition on the wetted surface and interior free surface as follows:

$$2\pi \sigma (P)+{\displaystyle \underset{H\cup F\prime }{\iint }\sigma (Q)G_n (P,Q) \mathrm{d}S}=\left\{\begin{array}{ll}{n}_j& j=1\dots 6\\ -\partial {\varphi }_0/\partial n& j=7\end{array}\right.,$$

(6)

$$4\pi \sigma (P)-{\displaystyle \underset{H\cup F\prime }{\iint }\sigma (Q)G_n (P,Q) \mathrm{d}S}=0.$$

(7)

It should be noted that ${\varphi }_n^{\prime }$ presents the normal derivative of the velocity potential in the interior computational domain. For points on the wetted surface located at a singular point on the hull ($P=Q$), the gradient of the Green function becomes singular and the small area around the singular point is excluded from the integration. The first term of the left-hand side of Equation (6) contributes to the value of normal derivatives in the vicinity of the singular point on the hull. The reference coordinate system for the mesh, as well as ship coordinate system, are presented in Figure 1.

As previously mentioned, added resistance in waves is the time-averaged second-order force, whose constant part can be determined based on the first-order potential. This can be achieved by integrating the first-order pressure over the mean wetted surface and considering the variation in the first-order quantities caused by the ship motion. As a response to the incoming wave, a ship is shifted from its equilibrium position resulting in the modified pressure distribution on wetted surface. Represented by the diagonal members of quadratic transfer function only, the added resistance for certain incoming wave frequencies can be determined based on ship response, the first-order velocity potential, the gradient of the velocity potential, the wave elevation on mean wetted surface, and along the waterline as follows [23]:

$$F=-\frac{1}{2}\rho g{\displaystyle \underset{WL}{\oint }{\zeta }_r^2}n \mathrm{d}l+\frac{1}{2}\rho {\displaystyle \underset{H}{\iint }\nabla \mathrm{\Phi }\nabla \mathrm{\Phi }} n \mathrm{d}s+\rho {\displaystyle \underset{H}{\iint }X} \cdot \nabla {\mathrm{\Phi }}_t n \mathrm{d}s+m\cdot R\cdot {\ddot{X}}_G,$$

(8)

where ${\zeta }_r^{}$ is the relative wave elevation, and $X=X_G+R\cdot x$ is first order motion vector in the ship coordinate system, determined based on the position of vector $x$ and the linearized rotation transformation matrix $R$, while ${\ddot{X}}_G$ is the acceleration vector in the ship’s centre of gravity. The angle of hull emergence is accounted for by modifying the normal vector as $n\prime =n/\left|\cos \gamma \right|$. The relative wave elevation is calculated based on the wave elevation along the hull, which is determined by the pressure on the panels along the waterline using the kinematic-dynamic boundary condition on the free surface.

As larger ships are being built, the wavelength to ship length ratio is becoming smaller, emphasizing the importance of estimating the added resistance in short waves as accurately as possible. In order to improve the numerical results obtained using the 3D panel method in short waves, the correction according to [25] is applied:

$$R_{AWr}=\frac{1}{2}\rho g{\zeta }_a^{2}B{B}_f{\alpha }_d(1+{\alpha }_U),$$

(9)

where $B_f$ is the bluntness coefficient obtained by the integration along the non-shaded part of the waterline, ${\alpha }_d$ is a term based on the draught and encountering wave number, and ${\alpha }_U$ is a term based on the forward speed. Finally, the added resistance in regular waves is defined as:

$$R_{AW}=R_{AW m}+R_{AW r}.$$

(10)

The coefficient of added resistance in waves is obtained as follows:

$$C_{AW}=\frac{R_{AW}}{\rho g{\zeta }_a^2{B}^2/L},$$

(11)

where $B$ is the ship breadth, $L$ the ship length, and ${\zeta }_a^{}$ is the wave amplitude. The added resistance in irregular waves is calculated on the basis of the spectral analysis. Namely, using the encounter wave energy spectrum $S_{\zeta }\left({\omega }_e\right)$ and added resistance in regular waves, the averaged value of added resistance in a frequency range is calculated by:

$$\overline{R_{AW}}=2{\displaystyle \underset{0}{\overset{\infty }{\int }}{S}_{\zeta }\left({\omega }_e\right)\frac{R_{AW}}{{\zeta }_a^2}\left({\omega }_e\right)\mathrm{d}{\omega }_e},$$

(12)

where ${\omega }_e=\omega -\frac{{\omega }^2}{g}v\cos \beta$ is the encounter frequency dependant on the wave heading $\beta$. Within this research, the two-parameter Bretschneider wave energy spectrum, as the modification of the one-parameter Pierson Moskowitz spectrum, was used to represent fully developed seas, while the Joint North Sea Wave Project (JONSWAP) was used to represent long crested limited fetch seas, as recommended by International Towing Tank Conference (ITTC) [26]. The added resistance at different sea states using the Bretschneider and JONSWAP sea spectra is calculated as follows:

$$S_{\zeta B}\left(\omega \right)=\frac{173H_{S1/3}^2/{\left(1.086 {\overline{T}}_Z\right)}^4}{{\omega }^5}{\mathrm{e}}^{-\frac{692/{\left(1.086 {\overline{T}}_Z\right)}^4}{{\omega }^4}},$$

(13)

$$S_{\zeta J}\left(\omega \right)=\frac{320H_{S1/3}^2}{T_p^4{\omega }^5}{\mathrm{e}}^{\frac{-1950}{T_p^4{\omega }^4}}{3.3}^{e^{\left\{-{\left(\frac{\omega /{\omega }_p-1}{\sigma \sqrt{2}}\right)}^2\right\}}},\hspace{1em}\hspace{1em}\hspace{1em}\sigma =\left\{\begin{array}{l}0.07\hspace{1em}\mathrm{for}\hspace{1em}\omega <{\omega }_p\\ 0.09\hspace{1em}\mathrm{for}\hspace{1em}\omega >{\omega }_p\end{array}\right.,$$

(14)

where $H_{S1/3}^{}$$\left(H_S\right)$ is the significant wave height, ${\overline{T}}_Z$ is the zero crossing period, and $T_p^{}=\frac{2\pi }{{\omega }_p}$ is the peak period defined with respect to the zero crossing period as $1.073{\overline{T}}_Z=0.834T_p$.

Validation and Verification Studies

The numerical results obtained using the 3D panel method were validated against the experimental data for the KRISO container ship (KCS) at 24 knots, corresponding to Froude number $Fn$ = 0.26. The main characteristics of the KCS are shown in

Table 1. The main particulars of KCS.

Main Particular	Value
$L_{PP}$, m	230
$B_{WL}$, m	32.2
$T$, m	10.8
$\nabla$, m3	52,030
$C_B$	0.651
$S$, m2	9530
$LCB$, %	−1.48
$LCG$, m	111.6
$VCG$, m	7.28
$r_{xx}/B$	0.40
$r_{yy}/L$	0.25
$V$, kn	24

. Relative deviation between the numerical and experimental results is calculated as follows:

$$RD=\frac{C_{AW,\mathrm{NUM}}-C_{AW,\mathrm{EXP}}}{C_{AW,\mathrm{EXP}}}100\%.$$

(15)

The verification procedure for the KCS was carried out according to [27] to assess the influence of the number of ship hull panels on the added resistance in waves, and to quantify the numerical uncertainty of the obtained results. The grid refinement was performed by subdividing each panel into four panels, resulting in three panel models with 1048, 4366, and 17,823 panels on the hull. The refinement coefficient for the panels on the interior free surface was equal to 1.5 and was kept constant. Coarse and fine panel models can be seen in Figure 2.

The assessment of the numerical uncertainty of the added resistance coefficient in regular and irregular head waves was carried out according to procedure recommended by the ITTC [28] and applied in [29]. The numerical uncertainty of added resistance in irregular waves was assessed for sea state characterized by $H_s=3.5 \mathrm{m}$ and $T_Z=10.5 \mathrm{s}$ for Bretschneider wave energy spectrum.

The type of convergence can be determined as follows [28]:

$$R=\frac{\left(S_2-S_1\right)}{\left(S_3-S_2\right)},$$

(16)

where $S_2-S_1$ is the difference between solutions obtained using a medium and fine grid, and $S_3-S_2$ is the difference between solutions obtained using a coarse and medium grid. For $0<R<1$, the monotonic convergence was obtained, and $-1<R<0$ denotes oscillatory convergence, while $\left|R\right|>1$ denotes divergence. Using the generalized Richardson extrapolation (RE), numerical uncertainty can be determined in the case of monotonic convergence. The order of accuracy reads:

$$p=\frac{\ln \left(\left(S_3-S_2\right)/\left(S_2-S_1\right)\right)}{\ln \left(r\right)},$$

(17)

where $r$ is the refinement ratio calculated using the total number of panels $h_i=1/\sqrt{N_i}$ [30,31]. RE error is defined as:

$${\delta }_{\mathrm{RE}}=\frac{S_2-S_1}{{\left(\frac{h_2}{h_1}\right)}^p-1}.$$

(18)

The normalized numerical uncertainty can be determined using the safety factor $F_S=1.25$ as follows:

$$\overline{U}=\frac{F_S \left|{\delta }_{\mathrm{RE}}\right|}{S_1-{\delta }_{\mathrm{RE}}}\cdot 100$$

(19)

3. Data Acquisition

In this study, the possibility of the determination of added resistance is investigated for container ships. Container ships were selected because, as one the largest seaborne vessels transporting over 90% of non-bulk cargo, they emit about 23% of the total CO₂ emissions caused by international maritime transport [32]. Container ships operate on established sailing routes between the busiest and the best connected container ports. Due to the scheduled time of arrival in ports, they typically sail at higher speeds compared to other types of merchant ships. Moreover, they are independent of market demands due to the variety of cargo they transport, and, for this reason, the number and size of container ships is continuously growing. For example, in the period 2015–2018 the number of container ships increased by approximately 5.1%, while the number of ships in the world fleet increased by approximately 8.5% [33], causing an increase in CO₂ emission of about 6.5% [34]. Hull forms of modern container ships with different dimensions, bow and stern types, section type, and block coefficient (prismatic coefficient) were generated in order to create a sufficiently large database for ANN training. The original 15 hull forms [35] were modified, resulting in an additional 75 hull forms with different prismatic coefficients and longitudinal centre of buoyancy positions (

Table 2. The original hull forms of container ships.

Ship	${\mathit{L}}_{\mathit{P}\mathit{P}}, \mathbf{m}$	$\mathit{B}$, m	$\mathit{T}$, m	$\mathit{L}\mathit{C}\mathit{B},\%$	${\mathit{C}}_{\mathit{M}}$	${\mathit{C}}_{\mathit{B}}$	${\mathit{C}}_{\mathit{P}}$	$\nabla$, m3
CS1	319	42.8	13	−3.49	0.9682	0.5907	0.6100	104,841.9
CS2	153	25	8	−0.691	0.9836	0.6583	0.6693	21,357.9
CS3	132	21.5	7	−0.691	0.9836	0.6583	0.6693	13,865.9
CS4	247	32.26	12	−2.438	0.9589	0.6057	0.6316	61,018.5
CS5	178	25.85	9	−0.849	0.9598	0.5853	0.610	24,457.2
CS6	234	32.42	11.25	−1.464	0.9506	0.6812	0.7166	59,445.3
CS7	106	20.3	4.25	0.138	0.9622	0.7146	0.7426	6221.82
CS8	104.8	18	7.9	−0.499	0.9697	0.6173	0.6366	9937.3
CS9	155.4	23.3	9.2	−1.683	0.9878	0.5619	0.5689	19,473.8
CS10	300	37	11	−2.963	0.9305	0.5292	0.5687	67,244
CS11	196	28	8.5	−1.614	0.9622	0.5721	0.5946	26,687.6
CS12	120.7	19	7.5	−0.886	0.9824	0.6682	0.6802	12,537.6
CS13	215	27	8	−0.216	0.9769	0.7097	0.7265	34,808
CS14	360	49	15	−1.137	0.9632	0.5830	0.6053	164,760.5
CS15	340	44	15	−2.438	0.9588	0.6100	0.6360	143,199.6

). The input variables for ANN were selected based on the analysis of the influence of ship characteristics on added resistance performed within [13]. It should be noted that the input variables should be determined in such a way that their influence on the output is more significant than their numerical uncertainty, bearing in mind that the input vector should be parsimonious. The influence of the ship mass characteristics, i.e., the vertical position of the centre of gravity and the pitch radius of gyration, on the added resistance in regular waves was investigated in [13,36], while the ship design speed was determined based on the results of regression analysis available in the literature [37,38]. The results of the regression analysis also defined the variation range of the prismatic coefficient with respect to the limitations of the method used for the hull form modification [39]. By changing the longitudinal position of the ship cross sections, keeping the same shape and area of cross sections, and the original midship coefficient, fore and aft prismatic coefficients were changed, resulting in either the change in the total prismatic coefficient or the longitudinal position of the centre of buoyancy.

It should be noted that hydrodynamic calculations of the added resistance in regular waves using the 3D panel method were performed for pitch radius of gyration equal to 24%, 25%, and 26% of $L_{PP}$, while the roll radius of gyration was equal to 35% of $B$. By performing the numerical calculations, 65,736 samples in total were obtained: 70% of samples (46,015) were used for training, 25% (16,434) for validation, and 5% (3,287) for testing ANN. During the training and validation processes, the ANN was familiar with both input and output values, while for the testing purposes the ANN had only input values at its disposal. The ranges of input variables were the following: the length between perpendiculars $L_{PP}$ = 104.8 ÷ 360 m, breadth $B$ = 18 ÷ 49 m, draught $T$ = 4.25 ÷ 15, displacement volume $\nabla$ = 5648.7 ÷ 174,090 m³, longitudinal position of the centre of buoyancy $LCB$ = 49 ÷ 180.72 m, block coefficient $C_B$ = 0.510 ÷ 0.780, prismatic coefficient $C_P$ = 0.530 ÷ 0.811, speed $V$ = 12.1 ÷ 26 knots, and pitch radius of gyration $r_{yy}$ = 25.15 ÷ 93.6 m. The top container ship sailing routes were analysed and the sea states with the highest probability of occurrence for waves coming from all directions and head waves for a given route were extracted based on the global wave atlas. In this way, the results of added resistance in irregular waves i.e., at different sea states that ship may encounter during sailing, were obtained. One of the input variables of the ANN is the choice of wave energy spectrum. Ranges of significant wave height and zero crossing period as input variables were $H_S$ = 0.5 ÷ 8.5 m, and $T_Z$ = 3.5 ÷ 12.5 s, respectively. It should be noted that some of the obtained numerical results were discarded if $H_S$ was larger than $T$.

4. Artificial Neural Network Model

In this section, a mathematical model of the feedforward MLP neural network, with one hidden layer and an error backpropagation (EBP) algorithm used in the estimation of the added resistance for actual sea states, is presented. The reasons for this was the possibility of supervisory training and relatively simple implementation. Based on the analysis of different neural network topologies, i.e., by varying the number of neurons in the hidden layer, learning rate and momentum parameters, as well as by performing an analysis of the results obtained using different learning algorithms, an adequate artificial neural network was proposed. The training process was performed by using the results of the hydrodynamic calculations, which were adequately prepared and divided into data sets for training, validation, and testing. The learning and testing accuracy of the neural network were evaluated based on the normalized value of the root mean square error (NRMSE), which is calculated as follows:

$$\mathrm{NRMSE}=\frac{\mathrm{RMSE}}{{\sigma }_{d_n}}=\frac{\sqrt{\frac{{\displaystyle \sum _{n=1}^N{\left(d_n-O_n\right)}^2}}{N}}}{{\sigma }_{d_n}},$$

(20)

where ${\sigma }_{d_n}$ is the standard deviation, $d_n$ is the desired value of output, $O_n$ is the output given by the ANN, and $N$ is the number of learning, testing, and validation samples. An example of the MLP neural network with a single hidden layer is given in Figure 3. The values of the input variables were multiplied by parameters (or weights), and an activation function in the hidden neurons was applied to a weighted sum. A bias neuron of the constant unit value was added to both the input and hidden layers, providing an additional degree of freedom in the learning process of ANN. The activation function in hidden neurons was selected as a unipolar sigmoid, which was a differentiable and bounded function with large values of derivatives near the origin and saturation towards infinity:

$$\gamma =\frac{1}{1+e^{-x}}.$$

(21)

The activation function of the output neuron was set as linear, allowing the ANN to learn nonlinear and linear relationships between input and output vectors, and giving values outside the range of −1 to +1. In order for the ANN to approximate the given nonlinear regression function, optimal weights should be applied. The aim of the training process was to iteratively determine the learning parameters and weights by minimizing the performance function based on one of the learning algorithms. The backpropagation algorithm used the samples of the inputs and corresponding desired outputs and calculated the errors in the forward direction. The initial weights were set as random values of the normal distribution in a range $\pm w_{\max }=\sqrt{3/\left(n-1\right)}$, where $\left(n-1\right)$ corresponded to the number of input variables, and the weights were iteratively adjusted to minimize the error function. When the steepest gradient descent (SGD) learning algorithm was used, the weights $\vartheta$ were changed based on the negative value of the learning error gradient of the performance function, as follows [40]:

$$\vartheta \left(n+1\right)=\vartheta \left(n\right)-\eta \nabla E\left(\vartheta \left(n\right)\right)+\alpha \Delta \vartheta \left(n-1\right),$$

(22)

where $n$ is the iteration step, $\eta$ is the learning rate, $\nabla E$ is a gradient of the performance function, $\alpha$ is the first order momentum, and $\Delta \vartheta$ is the change in weights in the previous iteration step. The performance function $E$ was calculated as the sum of the squared difference between the desired output values and output values obtained by the ANN for all training samples:

$$E=\frac{1}{2}{\displaystyle \sum _{n=1}^N{\left(d_n-O_n\right)}^2}.$$

(23)

To improve the convergence rate of the training or learning procedure, in this paper the learning rate was used and set as a constant value equal to $\eta =0.01$, based on the evaluated error of the training and validation data sets (Figure 4). It was evident that the slope of the error function decreased and the error of the validation data set increased for $\eta >0.01$. The ANN used to determine the optimal value of the learning rate had 55 neurons in the hidden layer, and this number was derived by analysing the NRMSE of training (TNRMSE) and validation (VNRMSE) data sets (Figure 5). Namely, the ANN topology with 55 neurons in the hidden layer, TNRMSE and VNRMSE, had the lowest values, which were obtained without the momentum included. Momentum, used to accelerate the training, is proportional to the weight update at the previous iteration step. It enables training in flat regions of the error surface and avoidance of the local optima, but sometimes it does not ensure learning error convergence. While weights are optimized through an iterative process of minimizing the error function, the ANN topology was determined experimentally. The possibility of applying the variable learning rate was also investigated. Namely, based on the comparison between the error in current and previous iteration steps, if the error in the current iteration step was larger than the error from the previous step (increased by 4%), the learning rate was decreased for 10%, and the weight update was discarded. Otherwise, the weights were updated, and the learning rate was increased by 1%. It should be emphasized that the mentioned rates of increase or decrease of the learning rate were determined experimentally.

In multidimensional nonlinear mapping problems, the error surface can be extremely complex with numerous local optima and flat regions, and the curvature may not be the same in all directions. Thus, the error surface can have a large slope in one direction, not necessarily leading to minimum point, and a small slope in the other direction, with the gradient pointing in that direction being very small. Therefore, it can be beneficial to use the information about the curvature along with the gradient of the error surface. For that reason, learning algorithms based on the second order derivatives were employed as well. However, since the exact calculation of the second order derivatives was computationally quite expensive, the learning algorithms that approximated the second order derivatives were employed. One of the second order learning algorithms analysed within this paper was the Levenberg-Marquardt (LM) algorithm [41], which corresponds to the gradient descent method, when the weights are far from their optimal values, and to the Gauss-Newton method, when the weights are getting closer to their optimal values. Namely, the Gauss-Newton method assumes the error function locally quadratic and finds the minimum. The Hessian matrix was approximated using the Jacobian matrix $J\left(n\right)$ that contains errors first derivatives with respect to the weights. The LM algorithm updated the weights according to the following equation:

$$\vartheta \left(n+1\right)=\vartheta \left(n\right)-{\left[J^T\left(n\right)J\left(n\right)+\mu I\right]}^{-1}{J}^T\left(n\right)e\left(n\right)$$

(24)

As a variation of a conjugate gradient method, the scaled conjugate gradient (SCG) algorithm [42] was applied as well. Unlike other conjugate gradient methods, the SCG algorithm does not require the line-search technique in each iteration step, but uses a Levenberg-Marquardt approach for the scale of the step size. In the first iteration step, the algorithm calculates the negative gradient $G=-\nabla E$ to start the search. The search direction in each next iteration step was conjugated to the previous search direction, i.e., the weights were updated according to the conjugate direction $D$ as follows:

$$\vartheta \left(n+1\right)=\vartheta \left(n\right)+\eta (n)D(n)$$

(25)

The conjugate direction at each iteration step was calculated as follows:

$$D\left(n+1\right)=G\left(n+1\right)+\frac{\left({\left|G\left(n+1\right)\right|}^2-G{\left(n+1\right)}^{T}G\left(n\right)\right)}{D{\left(n\right)}^{T}G\left(n\right)}D\left(n\right)$$

(26)

To improve the generalizability of the developed ANN model, the Bayesian regularization (BR) [43] was applied as well: the ANN with the BR is more robust than the standard EBP neural network, without the need for validation, which is one of its greatest advantages. The BR converts a nonlinear regression problem into a statistical problem with a probabilistic approach. They are difficult to overfit, since the training is performed with a limited number of effective weights, discarding the ones that are not relevant. The optimal weights, as well as the effective number of weights, are obtained by maximizing the posterior probability. Since the BR requires the determination of the Hessian matrix, it was approximated using the Gauss-Newton method within the LM learning algorithm. Alongside the sum of the squared error of the predicted and desired outputs, the error function included the sum of squared weights as follows:

$$F=\beta E_D+\alpha E_W,$$

(27)

where $\beta$ and $\alpha$ are the so-called regularization parameters that can be determined as follows:

$$\alpha =\frac{N-2\alpha \mathrm{tr}{\left[{\nabla }^{2}F\left(\vartheta \right)\right]}^{-1}}{2E_W\left(\vartheta \right)},$$

(28)

$$\beta =\frac{n-\left(N-2\alpha \mathrm{tr}{\left[{\nabla }^{2}F\left(\vartheta \right)\right]}^{-1}\right)}{2E_D\left(\vartheta \right)},$$

(29)

where $n$ is the number of training samples, $N$ is the number of weights, and the Hessian matrix is approximated as ${\nabla }^{2}F\left(\vartheta \right)\approx 2\beta J^{T}J+2\alpha I_N$.

More details regarding the mathematical background of the applied learning algorithms can be found in the before cited literature.

Preparation of Data

To ensure that input variables of different ranges and dimensions had an equivalent impact on output variables, the data in the learning set needed to be standardized. Otherwise, the input variables of the smallest values would have a negligible impact on the output variables and would not be accounted for during the learning process. The standardization of the input variables ensures a unit standard deviation and the mean value of the particular input variable data equal to zero. Even though the input data do not follow the normal distribution, they are standardized according to following equation:

$$Z_i=\frac{Z_i{}^{\prime }-\frac{1}{N}{\displaystyle \sum _{i=1}^N{Z}_i{}^{\prime }}}{\sigma }.$$

(30)

where $\sigma$ is a standard deviation. For comparison purposes, the input data was also normalized prior to training, i.e., the range of data was between 0 and 1. It should be mentioned that the output data were also normalized.

To accelerate training and investigate its influence on the residuals, principal component analysis (PCA) was performed for the standardized input data, which showed different dynamics of the change in the data. Namely, the PCA was a dimensionality reduction method, applied to remove the collinearity of data prior to training, thus enabling the reduction in the number of input variables, while preserving as much data variance as possible (i.e., the information contained in the initial variables). Within this research, the number of new variables, i.e., principal components ${\xi }_p$, corresponded to the number of initial variables, but they were linearly independent, and are obtained by eigenvalue decomposition of the data covariance matrix as follows:

$${\xi }_p=w_{p1}{x}_1+w_{p2}{x}_2+\cdots +w_{pp}{x}_p,$$

(31)

where $x_p$ corresponds to the input variables, and $w_{pp}$ to the eigenvector of the data covariance matrix. It should be noted that the order of eigenvectors in the matrix is based on their eigenvalues.

When calculating the added resistance in waves for different sea states, different noise level in the data was observed depending on the frequency range of wave energy. For better fitting, the data were classified into three classes based on the wave period, which significantly improved the accuracy of the model as will be explained later.

5. Results and Discussion

As previously mentioned, the hulls of container ships are discretized using quadrilateral panels. Each ship is discretised by approximately 4000 panels on the wetted surface and 6000 panels on the interior free surface, with the refinement coefficient equal to 0.5 for the irregular frequency treatment. Numerical calculations were performed for full-scale ships, in the frequency range of incoming regular head waves from 0.1 to 1.5 rad/s with a step of 0.05 rad/s. The added resistance was determined based on the direct pressure integration in unrestricted water.

5.1. Validation and Verification of the Numerical Results

In comparison to the available experimental data, the numerically obtained coefficient of the added resistance in waves for KCS ships, with the correction in short waves applied, showed satisfactory agreement, with the maximum relative deviation (RD) equal to 8% (Figure 6). The numerically obtained results significantly underestimated the added resistance in short waves. Waves are considered short relative to ship length when the ratio between the wavelength and ship length is $\lambda /L<0.5$. However, the influence of wave diffraction in short waves can be observed for $\lambda /L>0.5$ as well, and the obtained results, without the correction, showed large relative deviation compared to the experimental results. Similar occurrence was observed in the case of other analysed hull forms of container ships. As can be seen from Figure 6, large discrepancies between the experimental results in the high frequency region were present as well.

The results of the performed verification study for panel size used to discretise KCS ships are shown in

Table 3. Results of the verification study for panel size used to discretize KCS ships.

$\mathit{\lambda }/\mathit{L}$	${\mathit{S}}_3$	${\mathit{S}}_2$	${\mathit{S}}_1$	$\mathit{R}$	$\mathit{p}$	${\mathit{\delta }}_{\mathit{R}\mathit{E}}$	${\mathit{S}}_0$	$\overline{\mathit{U}},\%$
2.0	1.518	1.513	1.492	4.4882	/	/	/	/
1.5	6.326	6.244	6.221	0.2911	1.7803	0.0098	6.211	0.197
1.33	9.193	9.069	8.989	0.6483	0.6252	0.1484	8.840	2.098
1.15	11.279	10.996	10.718	0.9828	/	/	/	/
0.5	0.780	0.762	1.568	−43.1464	/	/	/	/
$C_{AW}$	10.813	10.578	10.455	0.5238	0.9329	0.1353	10.320	1.639

. The monotonic convergence, in the case of regular waves, was obtained for only $\lambda /L=1.5$ and $\lambda /L=1.33$. The numerical uncertainty of the obtained added resistance coefficient in regular waves was lower than 2.1%. The coefficient of the mean added resistance in irregular waves can be seen in Table 3. The monotonic convergence was obtained and the calculated numerical uncertainty amounted to 1.64%. It should be noted that the number of panels for coarse, medium and fine panel models equals 1048, 4366 and 17,823, respectively.

5.2. Results of Artificial Neural Network

As previously mentioned, the input variables of the ANN model were determined based on the results of the sensitivity analysis presented in [13]. The ANN model predicted the value of the added resistance based on 12 input variables: $R_{AW}=f(L_{PP},B,T,\nabla ,LCB,C_B,C_P,V,r_{yy},S_{\zeta },T_Z,H_S)$. The PCA was performed for the standardized input data, with the mean value equal to zero and a standard deviation equal to one. Based on the covariance matrix, the linear dependence between the input variables can be observed, except in the case of wave energy spectrum $S_{\zeta }$, which was equal to 1 for Bretschneider and to −1 for the JONSWAP spectrum, and $T_Z$ and $H_S$ with the ship characteristics and speed. On the other hand, the relation between $T_Z$ and $H_S$ was noted. The percentage variance of each principal component was calculated as the ratio between the variance of the principal component and the total variance (Figure 7).

It can be seen that most of information was contained in the first principal component (PC), and, that by eliminating the last five PCs with the smallest variance, 99.4% of the total variance was retained. However, as already mentioned, within this study, all PCs were accounted for in the learning process, as the aim of the PCA was not dimensionality reduction. A significantly shorter computational time was observed when the ANN was trained with the PCs. When comparing TRNMSE and VRNMSE during the training with the initial input variables and the PCs, using the gradient descend learning algorithm, a larger decrease in the error after the same number of iteration steps was observed in the case of the PCs.

The proposed ANN model, with the error back propagation learning procedure, was trained with normalized and standardized input data, and PCA data as well. The selected initial weights were kept the same for the initialization of the learning process using different learning algorithms. The coefficients obtained within the PCA were used to prepare both the data for training and validation. Different momentum values within the SGD algorithm were analysed for the constant value of the learning rate equal to 0.01. The obtained results are presented in

Table 4. The obtained results for training and validation errors for different learning algorithms.

Learning Algorithm	SGD		LM		SCG		BR
Error	TNRMSE	VNRMSE	TNRMSE	VNRMSE	TNRMSE	VNRMSE	TNRMSE
Norm. data 1	0.0624 $(\alpha =0$)	0.0766 $(\alpha =0$)	0.0231	0.0239	0.0583	0.0589	0.0231
	0.0615 $(\alpha =0.3$)	0.0756 $(\alpha =0.3$)
	0.0613 $(\alpha =0.5$)	0.0736 $(\alpha =0.5$)
	0.0647 $(\alpha =0.8$)	0.0703 $(\alpha =0.8$)
Stand. data 2	0.0633 $(\alpha =0$)	0.0772 $(\alpha =0$)	0.0238	0.0246	0.0581	0.0585	0.0238
	0.0616 $(\alpha =0.3$)	0.0757 $(\alpha =0.3$)
	0.0618 $(\alpha =0.5$)	0.0742 $(\alpha =0.5$)
	0.0665 $(\alpha =0.8$)	0.0725 $(\alpha =0.8$)
PCA data	0.0639 $(\alpha =0$)	0.0792 $(\alpha =0$)	0.0234	0.0236	0.0545	0.0559	0.0204
	0.0628 $(\alpha =0.3$)	0.0788 $(\alpha =0.3$)
	0.0619 $(\alpha =0.5$)	0.0774 $(\alpha =0.5$)
	0.0670 $(\alpha =0.8$)	0.0804 $(\alpha =0.8$)

¹ Normalized data. ² Standardized data.

. It should be noted that TNRMSE and VNRMSE for the SGD algorithm were obtained after 10,000 iteration steps, for the SCG algorithm after 6578 iteration steps, for the LM algorithm after 1884 iteration steps, and for the LM algorithm with BR after 5822 iteration steps. It is important to mention that the training was stopped when the validation error began to stagnate. It was observable that there were no significant differences between the results obtained after training the ANN with normalized and standardized input data. In the case of the SGD and LM learning algorithms, slightly lower values of TNRMSE and VNRMSE were obtained when trained with normalized data. The ANN trained with the PCs indicated the preservation of the accuracy, while reducing the computational time significantly, despite the fact that all the PCs were kept. The PC are linear combinations of the initial variables, but completely linearly independent of each other. Based on the obtained results using the SGD algorithm, it can be seen that the error of the validation data set was of the same order of magnitude as the error of the learning data set, and the change in momentum did not contribute significantly to the acceleration of the learning process or the improvement of accuracy. Namely, when the momentum was applied, the weights updated depending on the error gradient from the previous iteration steps, and as the number of steps increased, the influence of the previous weights update contained in the momentum term was reduced. In complex multidimensional nonlinear problems, the error surface at a certain point can have a small slope in one direction, although it is the right direction towards minimum, and large slope in the other direction, leading to a flat area or local minimum of the error surface. Since the weights are updated by the largest negative error gradient, the search direction towards the global optimum changes. Therefore, in addition to the error gradient of the first order, the results obtained using learning algorithms that approximate the second order derivatives, were analysed as well, and were significantly improved compared to the SGD (Table 4). A variable learning coefficient is generally considered more effective compared to using a constant value for the learning rate during training. However, using the variable learning rate can decrease the training error, but increase the validation error leading to poor the generalizability of the ANN. The application of a variable learning rate led to a decrease of the TNRMSE (0.0476), but also a significant increase in the VNRMSE (0.1568).

There are no large deviations between TNRMSE and VNRMSE obtained using the LM and SCG algorithms, and the LM algorithm with Bayesian regularization (BR) trained with normalized, standardized or PCA data, although the error seemed to be the smallest in the case of training with PCA data. The lowest computational time was achieved for the SCG algorithm, but the accuracy of that ANN was also the lowest. The ANN with the LM algorithm trained with PCA data showed slightly better generalizability compared to the results obtained using the normalized and standardized data.

The ANN trained by the LM algorithm with BR proved to be the most effective in the estimation of the added resistance in waves. Although the required computational time was slightly higher compared to the other algorithms, with the exception of the SGD, the TNRMSE was equal to 0.0204 in the case of training with PCA data, which was the lowest error obtained. Since BR does not require the validation set, all data were used for training. In other words, BR allowed 95% of the data to be used for the training. The value of the learning rate for the LM algorithm was set to 0.004, and 751 out of a total of 771 weights were effectively involved in the learning process. It can be concluded that the selected network topology with 55 neurons in the hidden layer was adequate to describe the given hydrodynamic problem. Since the number of effective weights was close to the total number of weights, 65 neurons in the hidden layer were employed. The TNRMSE obtained with this ANN topology was equal to 0.0183 after 2450 iteration steps, and 890 out of a total of 991 weights effectively participated in the learning process.

The regression plots of the results obtained with the SGD, SCG, LM and LM with BR algorithms, along with the coefficients of determination, for the network with 55 neurons in the hidden layer, are presented in Figure 8. Some deviations can be observed in the case of the SGD and SCG algorithms, especially for smaller values of added resistance in waves, while the LM with and without BR allowed for a more accurate estimate. It should be noted that all presented results were obtained using PCA data for training. Although the coefficients of determination were relatively high and the TNRMS relatively low, it was necessary to analyse the residuals and evaluate the generalizability of the ANN based on the data that were not used within the training. For this reason, 5% of the total amount of data was used for the network testing. Figure 9 shows the residuals of the data used for training and testing of the ANN with the LM and BR, and 55 neurons in the hidden layer. The residuals were randomly distributed and centred around the zero value, indicating that the applied nonlinear regression model was adequate. The NRMSE of the data for testing was equal to 0.0227, while the RD was equal to −0.42%, and −0.1% for the data used for training and testing, respectively. Although the RD was very low, the residual analysis showed that the data were overestimated or underestimated, cancelling the RD and leading to a low mean RD for all values.

The regression plot and RD of the data used for the testing of the ANN with LM and BR with 55 neurons in the hidden layer are shown in Figure 10. The ANN with 65 neurons in the hidden layer had an RD equal to −0.46% and −2.33% for the data used for learning and testing, respectively. The NRMSE of the test data set was equal to 0.0201. The regression plot and RD of the data used for the testing of the ANN with LM and BR with 65 neurons in the hidden layer are shown in Figure 11. A significant deviation between the desired outputs and the ones given by the ANN were observed in the case of the lower values of added resistance in waves. The analysis of the obtained results showed that the largest RD were obtained for significant wave heights equal to 0.5 m. Certain deviations in the results also occurred at lower values of the wave period. Namely, small values of added resistance, when the wave energy was concentrated in the range of high frequencies, and when the frequency range of the response spectrum was not overlapping the frequency range of the wave energy spectrum, lead to small values of added resistance in certain sea states. Otherwise, the accuracy of the network was within 10%, meaning that, despite the different dynamics of change in the data and different noise levels depending on the frequency range of the response spectrum, the ANN can be considered successful in the estimation of added resistance in waves. The ANN model adapts very well to the nonlinearity of data, which can be seen for the first 100 samples of the testing data in Figure 12.

For analysed sea states, defined within

Table 5. The significant wave height and wave period for analysed sea states.

Sea State	${\mathit{H}}_{\mathit{S}}, \mathbf{m}$	${\mathit{T}}_{\mathit{Z}}, \mathbf{s}$
SS1	1.5	6.5
SS2	2.5	7.5
SS3	2.5	8.5
SS4	3.5	9.5
SS5	3.5	10.5
SS6	4.5	11.5

, absolute and percentage differences between the added resistance in waves obtained by hydrodynamic calculations and the ANN with the LM algorithm with BR, and 65 neurons in the hidden layer for the Bretschneider and JONSWAP wave energy spectra, are shown in Figure 13 and

Table 6. The absolute and percentage differences between the results obtained by ANN and BIEM for KCS.

Sea State	SS1	SS2	SS3	SS4	SS5	SS6
Bretschneider	−18,871.77 N (−29.74%)	−17,570.78 N (−9.32%)	−9839.05 N (−4.78%)	14,353.20 N (3.49%)	29,563.43 N (7.76%)	23,951.34 N (4.33%)
JONSWAP	4175.82 N (6.95%)	1610.69 N (1.00%)	−5053.55 N (−2.70%)	26,480.53 N (5.66%)	45,777.20 N (9.34%)	15,146.23 N (2.27%)

. It can be seen that the highest RD were obtained for SS1, while, for the other sea states, RD was below 10%. To develop a model based on the ANN that will adapt better to the data depending on the frequency range of the response spectrum, and to overcome the limited accuracy for lower values of added resistance, the data were divided into three classes based on the wave period. Namely, depending on the frequency range of the spectral energy for certain sea states, the influence of certain input variables on added resistance could be negligible or significant. Therefore, the first, second, and third classes cover the data with $T_Z=3.5÷6.5$ s, $T_Z=7.5÷8.5$ s, and $T_Z=9.5÷12.5$ s, respectively. The regression plot and RD of the data used for testing of ANN for three classes are shown in Figure 14. Again, larger RD could be observed for a significant wave height equal to 0.5 m and for smaller ships in short waves, for example a ship that is 132 m long for a sea state defined with a significant wave height of 3.5 m and a wave period of 6.5 s. Since the results of hydrodynamic calculations at high wave frequencies underestimated the added resistance due to the diffraction component, the results depended on the applied correction. Larger residuals do not necessarily mean that the accuracy of the model is limited, but may indicate good generalization properties of the network, since the ANN is not trying to adapt to the noise in the data. Similar observations can be made for the data in the second class as well. The TNRMSE was equal to 0.0112 and 0.0094 for the first and second classes, respectively. The NRMSE for the testing data was equal to 0.0116 and 0.0089 for the first and second class, respectively. The RD of the training and testing data were equal to −0.19% and 0.10% for the first class, and −0.76% and −1.36% for the second class, respectively. The TNRMSE of the third class was equal to 0.0080, and the NRMSE of the testing was equal to 0.0081. The RD was equal to 0.065% and −0.999% for the training and testing data, respectively. The ANN achieved satisfactory accuracy and generalizability, especially in the third class, with some more pronounced residuals corresponding to the smallest significant wave height, however, the absolute values of the added resistance in that case were relatively low. The limited accuracy of the ANN model in the estimation of added resistance in waves at lower sea conditions, when significant wave heights are around 0.5 m, or for very long waves, when spectral energy does not cover the response spectrum, was thus observed. This may be rationalized by the fact that, within the preliminary design, such small values would not affect the hull design or the propulsion system optimization.

Absolute and percentage differences between the added resistance in waves obtained by hydrodynamic calculations and the ANN, with the LM algorithm with BR and 65 neurons in the hidden layer for the Bretschneider and JONSWAP wave energy spectra after the data classification, are presented in Figure 15 and

Table 7. The absolute and percentage differences between the results obtained by ANN and BIEM for KCS after the data classification.

Sea State	SS1	SS2	SS3	SS4	SS5	SS6
Bretschneider	−979.25 N (−1.54%)	412.74 N (0.22%)	−10,368.66 N (−5.04%)	23,358.19 N (5.68%)	413.40 N (0.11%)	3841.04 N (0.69%)
JONSWAP	930.11 N (1.55%)	1258.65 N (0.78%)	8318.72 N (4.45%)	−1038.39 N (−0.22%)	−11,775.83 N (−2.40%)	−1692.20 N (−0.25%)

. It can be observed that the obtained results are of satisfactory accuracy, and that the numerical model based on the ANN has good generalizing properties when estimating the added resistance in waves for data that are not familiar to the ANN during training. Namely, as can be noticed from Table 7, the largest RD between the added resistance calculated numerically and the value obtained by the ANN was 5.68% for SS4 and the Bretschneider wave energy spectrum.

An ANN can be successfully used for the estimation of added resistance in waves since it enables fast estimation of ship added resistance in waves with satisfactory accuracy based on the ship characteristics, speed and sea state, which can contribute to an evaluation of added resistance of different hull forms in the design phase. ANN is a powerful tool, which, although not familiar with the physical background, can learn to describe the physical phenomenon and the relationship between input variables and the desired output. The developed numerical model has 12 input variables, 65 neurons in one hidden layer, and is based on the LM learning algorithm with BR. The model incorporated three sub models based on the classification of data according to the wave zero crossing period.

6. Conclusions

This paper investigated the possibility of using ANN for the estimation of ship added resistance in waves. The numerical model, based on the results of hydrodynamic calculations and ANN, can be used in the preliminary design phase of a ship with a block coefficient in the range of 0.510 and 0.780, or a prismatic coefficient in the range of 0.530 and 0.811, which sails at Froude numbers between 0.174 and 0.258. Hydrodynamic calculations were performed using BIEM in the frequency domain in regular head waves, and the results in short waves were corrected to account for the diffraction component adequately. The spectral analysis was performed using two theoretical wave energy spectra in order to determine the mean value of added resistance in actual sea states. The performed validation and verification studies for the KCS ship illustrated the satisfactory accuracy in comparison to the experimental data and low numerical uncertainty.

The accuracy and generalizability of ANN with different learning parameters and algorithms were investigated, and the LM algorithm with BR and 65 neurons in the hidden layer proved to be the most effective. Although a slightly more computationally demanding algorithm due to the determination and storage of Jacobi matrices, the LM learning algorithm with BR enabled the use of the most data for training, and a more robust model less sensitive to noise and outliers in the data was achieved. In addition to the advantages, the PCA performed using the standardized input data were pointed out, ensuring the same level of accuracy, but a significantly shorter computational time.

For different values of input variables, the results of the hydrodynamic calculations depended on the frequency range of the incoming regular waves, i.e., the frequency range of the wave energy spectrum. In that case, a complex physical model, which requires adequate preparation of a sufficient amount of data, is necessary. Thus, a way to potentially increase the accuracy and generalizability of ANN by data classification was presented within this paper. Such a developed numerical model allowed for an estimation of added resistance with satisfactory accuracy in a short time, depending on the class of sea state, without the need to perform complex hydrodynamic calculations, which could be of great practical use especially within preliminary ship design. Similar models could be developed for different ship types and various wave headings, which will form a part of future work, along with the analysis of advantages provided by other learning algorithms.