A Nonlinear Autoregressive Exogenous (NARX) Neural Network Model for the Prediction of the Daily Direct Solar Radiation

: The solar photovoltaic (PV) energy has an important place among the renewable energy sources. Therefore, several researchers have been interested by its modelling and its prediction, in order to improve the management of the electrical systems which include PV arrays. Among the existing techniques, artiﬁcial neural networks have proved their performance in the prediction of the solar radiation. However, the existing neural network models don’t satisfy the requirements of certain speciﬁc situations such as the one analyzed in this paper. The aim of this research work is to supply, with electricity, a race sailboat using exclusively renewable sources. The developed solution predicts the direct solar radiation on a horizontal surface. For that, a Nonlinear Autoregressive Exogenous (NARX) neural network is used. All the speciﬁc conditions of the sailboat operation are taken into account. The results show that the best prediction performance is obtained when the training phase of the neural network is performed periodically.


Introduction
The operation and management of the electrical system is becoming increasingly complex with the integration of renewable sources. The intermittent nature of renewable energies requires the prediction of the available power in the grid in order to maintain the production-consumption balance. This prediction is a complex task because it needs the forecasting of some variables and weather parameters, for instance, that of the direct solar radiation in the case of photovoltaic (PV) energy. Several methods for prediction of the solar radiation exist in the literature. They differ in the available inputs, their classification and the horizon of the prediction. Among these methods, some are based on linear models such as Linear Regression (LR), Auto-Regressive Moving Average (ARMA) and Auto-Regressive (AR) [1,2]. However, because of the nonlinear behavior of the solar radiation, researchers propose several nonlinear models based on wavelet-based methods, fuzzy models, Adaptive Neural Fuzzy Inference Systems (ANFIS) Random Forests (RF), k-Nearest Neighbors (kNN) and Artificial Neural Networks (ANN) [2][3][4][5].
Some researches combine linear and nonlinear methods such as in [6] where authors predict one minute ahead solar radiation using a hybrid method based on wavelet, ARMA and Nonlinear The necessity of the solar radiation prediction is related to several factors. On one hand, the solar radiation has an important influence on PV power production. On the other hand, it is greatly influenced by the sailboat mobility and the varying characteristics of the solar radiation during the year. In addition, it must be noted that this research study is a part of a bigger project which consists in providing electricity to a race sailboat using only renewable sources such as wind turbines, PV panels, hydro generator and energy recovery system. Thus, the nature of the project requires a high level of exactitude.
The solar radiation received at the Earth's surface depends not only on the Sun's position, but also on the optical state of the atmosphere. Therefore, the solar radiation model consists of two components: a deterministic component and a statistical component. The deterministic component depends only on the distance between the center of Sun and the point of measurement. It is easily identified by the laws of celestial mechanics and energetic physics as explained later in the manuscript [9][10][11][12][13][14][15][16][17][18][19]. In this research study, the deterministic component is calculated using the "clear sky" model. The statistical component depends on different factors: the cloud occurrence, the rate of ozone, the rate of humidity, etc. [20][21][22]. It is difficult to identify the statistical component due to the probabilistic feature of its different sub-components. Therefore, in this research study the statistical component is presented by its most important sub-component: the cloud cover.
Based on the above decomposition of the solar radiation, the "clear sky" model and the cloud cover are used as inputs for ANN to predict the daily direct solar radiation on a horizontal surface. The main contribution of this research study is the development of a new prediction model which takes into account two important factors: the change of the sailboat localization and the changing character of the solar radiation during the year such as the daily sunshine hours and the daily solar peak.
Considering the performance of the dynamic neural network in prediction issues and the advantages of the NARX models, the NARX-ANN is considered as prediction technique in this paper [4,23,24]. All neural network simulations are performed in a Matlab environment.
The rest of this paper is structured as follows: in Section 2, the technique used for the prediction is defined by describing the NARX network model. In Section 3, the model structure and the used database are introduced. Then, in Section 4 the deterministic component of the direct solar radiation is modelled. This section consists of two parts: the calculation of the geometrical parameters and the presentation of the direct solar radiation model. In Section 5, the statistical component of the direct solar radiation is presented. Section 6 focuses on the adjustments and the results and discussion. The last section concludes the paper.

Artificial Neural Networks and NARX Model
Artificial Neural Networks (ANNs) are used for an extensive range of issues as clustering, recognition, pattern classification, optimization, function approximation and prediction [25][26][27]. ANNs are mathematical tools stimulated to the biological neural system; they have indeed a powerful capacity to learn, store and recall information. They are black-box modeling tools which have the capacity to carry out non-linear mapping of an m-dimensional input space onto an n-dimensional output space, when the relations between the input and output spaces are unknown [28]. ANNs are developed based on mathematical models. The choice of the ANN model depends on the a-priori knowledge of the system to be modeled.
Considering that, on one hand, the solar radiation is a time series, and that, on the other hand, the NARX neural network is good predictor of time series [4,23,24], it is used in this research study. Actually, NARX concept is a nonlinear generalization of the Autoregressive Exogenous (ARX), which is a standard instrument in linear black-box system identification [29]. NARX models can be used to model a extensive variety of nonlinear dynamic systems. They have been applied in various applications including time-series modeling [30].
The NARX is a recurrent dynamic neural network. It has feedback connections which enclose several layers of the network. In order to obtain the full performances of the NARX neural network for Energies 2018, 11, 620 4 of 21 nonlinear time series prediction, it is interesting to utilize its memory ability using the past values of predicted or true time series.
As it can be seen in Figure 1, there are two different architectures of NARX neural network model, series-parallel architecture (named also open-loop) and parallel architecture (named also close-loop) given by the Equations (1) and (2), respectively: where, F(·) is the mapping function of the neural network,ŷ(t + 1) is the output of the NARX at the time t for the time t + 1 (it is the predicted value of y for the time t + 1).ŷ(t),ŷ(t − 1), . . . ,ŷ(t − n y ) are the past outputs of the NARX. y(t), y(t − 1), . . . , y(t − n y ) are the true past values of the time series, called also desired output values. x(t + 1), x(t), . . . , x(t − n y ) are the inputs of the NARX. n x is the number of input delays and n y is the number of output delays.
( ) ,..., 1 , , 1 , ,..., 1 , 1 (2) where, F(·) is the mapping function of the neural network, ŷ (t + 1) is the output of the NARX at the time t for the time t + 1 (it is the predicted value of y for the time t + 1). ŷ (t), ŷ (t − 1), …, ŷ (t − ny) are the past outputs of the NARX. y(t), y(t − 1), …, y(t − ny) are the true past values of the time series, called also desired output values. x(t + 1), x(t), …, x(t − ny) are the inputs of the NARX. nx is the number of input delays and ny is the number of output delays.
In the series-parallel architecture, the future value of the time series y(t − 1)is predicted from the present and past values of x(t) and the true past values of the time series y(t). In the parallel architecture the prediction is performed from the present and past values of x(t) and the past predicted values of the time series ŷ (t).
In this research study, during the training phase the series-parallel architecture is used because the availability of the true past values of the time series. The use of the series-parallel architecture has two advantages. The first is that the use of the true values as input of the feedforward network is more precise. The second advantage consists in an architecture of the resulting network which is purely feedforward, the usual training algorithms for Multi-Layer Perceptron (MLP) networks can be used. After the training phase, the NARX neural network is converted to the parallel architecture which is beneficial for multi-step-ahead prediction [31,32]. The mapping function F(·) is initially unknown and it is approximated during the training process of the prediction. In the NARX neural network model, the internal architecture that performs this approximation is the Multi-Layer Perceptron (MLP). The MLP offers a powerful structure allow to learn any type of continuous nonlinear mapping.
As shown in Figure 2, a classic MLP consists of three layers: input, hidden and output layer. Other elements consist of neurons, activation functions and weights. The direction of the information flow throughout the layers is from input to output layer. In each layer, each neuron multiplies the input vector xj given by the previous layer by the weights vector wij to give the scalar product × . An activation function f, given by Table 1, is then performed to obtain the following neuron output: In the series-parallel architecture, the future value of the time series y(t − 1) is predicted from the present and past values of x(t) and the true past values of the time series y(t). In the parallel architecture the prediction is performed from the present and past values of x(t) and the past predicted values of the time seriesŷ(t).
In this research study, during the training phase the series-parallel architecture is used because the availability of the true past values of the time series. The use of the series-parallel architecture has two advantages. The first is that the use of the true values as input of the feedforward network is more precise. The second advantage consists in an architecture of the resulting network which is purely feedforward, the usual training algorithms for Multi-Layer Perceptron (MLP) networks can be used. After the training phase, the NARX neural network is converted to the parallel architecture which is beneficial for multi-step-ahead prediction [31,32].
The mapping function F(·) is initially unknown and it is approximated during the training process of the prediction. In the NARX neural network model, the internal architecture that performs this approximation is the Multi-Layer Perceptron (MLP). The MLP offers a powerful structure allow to learn any type of continuous nonlinear mapping.
As shown in Figure 2, a classic MLP consists of three layers: input, hidden and output layer. Other elements consist of neurons, activation functions and weights. The direction of the information flow throughout the layers is from input to output layer. In each layer, each neuron multiplies the input vector x j given by the previous layer by the weights vector w ij to give the scalar product x j × w ij . An activation function f, given by Table 1, is then performed to obtain the following neuron output: where i is the index of the neuron in the layer. j represents the input index in the ANN. where i is the index of the neuron in the layer. j represents the input index in the ANN. The process of training consists in modifying the weights in an organized way using an appropriate algorithm. During the training process, a specified number of inputs and their desired output (named also target) are introduced in the network. Then, the weights are tuned so that the neural network produces an output close to the target values.
Another issue of the training phase is to deliver a globally optimal solution, avoiding local minima. For that, some weights are randomly initialized and the one which gives the best result is considered.

Model Structure and Used Database
As presented in Section 2, the prediction of a time series using ANNs requires two types of entries: the desired output of the ANN (named also the target) and a number of inputs which depends on the system to be modelled. Figure 3 shows the model structure defined in our study.
In our research study, the target is the measured global solar radiation on a horizontal surface. It is acquired from the weather station of "Ecole Supérieure des Technologies Industrielles Avancées" (ESTIA)-Bidart France. The measurements of the global solar radiation are performed with a step of 5 min time interval. It should be noted that measurements must be performed for the direct solar radiation. The weather station of ESTIA doesn't measure the diffuse solar radiation, so the calculation of the real direct solar radiation is not possible (Global radiation = Direct radiation + diffuse radiation). That's why the global solar radiation is taken as target of the NARX-ANN. Assuming that, the direct solar radiation consists of a deterministic component and a statistical component, the made model has two inputs:   The process of training consists in modifying the weights in an organized way using an appropriate algorithm. During the training process, a specified number of inputs and their desired output (named also target) are introduced in the network. Then, the weights are tuned so that the neural network produces an output close to the target values.
Another issue of the training phase is to deliver a globally optimal solution, avoiding local minima. For that, some weights are randomly initialized and the one which gives the best result is considered.

Model Structure and Used Database
As presented in Section 2, the prediction of a time series using ANNs requires two types of entries: the desired output of the ANN (named also the target) and a number of inputs which depends on the system to be modelled. Figure 3 shows the model structure defined in our study.
Energies 2018, 11, 620 6 of 21 calculated based on two coordinate systems: the equatorial coordinate system and the local coordinate system. The calculations are presented in Section 4.

•
The statistical component: is the exogenous input of the NARX model. Here, it contains only the cloud cover, because two reasons. First, the cloud cover is the most influential parameter on the direct solar radiation. Second, the identification of the other parameters is complicated. The construction of the cloud cover vector is presented in Section 5.

The Deterministic Component of the Direct Solar Radiation
This part presents the calculation of the deterministic component of the direct solar radiation.

Calculation of the Geometrical Parameters
Before presenting the direct solar radiation by a sophisticated model, it is appropriate to define the geometrical parameters describing the Sun position in a point of the Earth at a given moment. To reach this objective, two coordinate systems are used; the first is the equatorial coordinate system which defines the position of the Sun's center relative to the Earth's center. It is presented by the Sun's declination (δ) and the hour angle (ϖ). The second coordinate system is the local coordinate system which defines the position of the Sun's center relative to the point of measurement. It is presented by the sun height (α) and the azimuth angle (ψ).
This section focuses on the definition of all the geometric parameters required to calculate the direct solar radiation on a horizontal surface [33][34][35][36][37][38][39][40][41][42][43][44]. It should be noted that the literature contains more than one expression for some parameters. Therefore, a comparative study has been performed, in this section, in order to choose the adequate equations for three parameters: the Sun's declination (δ), the equation of time (E) and the Earth-Sun distance correction factor KD.

Sun's Declination
The Sun's declination (δ) is the first equatorial coordinate. It is the angle formed by the equator plane and the Earth-Sun direction ( Figure 4). The Sun's declination is approximated by Cooper as follows [33]: where, J is the rank of day in the current year (1 for January 1st  In our research study, the target is the measured global solar radiation on a horizontal surface. It is acquired from the weather station of "Ecole Supérieure des Technologies Industrielles Avancées" (ESTIA)-Bidart France. The measurements of the global solar radiation are performed with a step of 5 min time interval. It should be noted that measurements must be performed for the direct solar radiation. The weather station of ESTIA doesn't measure the diffuse solar radiation, so the calculation of the real direct solar radiation is not possible (Global radiation = Direct radiation + diffuse radiation). That's why the global solar radiation is taken as target of the NARX-ANN. Assuming that, the direct solar radiation consists of a deterministic component and a statistical component, the made model has two inputs: • The deterministic component: it is the endogenous input of the NARX model. It is described mathematically by the clear sky model of the direct solar radiation. The clear sky model is calculated based on two coordinate systems: the equatorial coordinate system and the local coordinate system. The calculations are presented in Section 4.

•
The statistical component: is the exogenous input of the NARX model. Here, it contains only the cloud cover, because two reasons. First, the cloud cover is the most influential parameter on the direct solar radiation. Second, the identification of the other parameters is complicated. The construction of the cloud cover vector is presented in Section 5.

The Deterministic Component of the Direct Solar Radiation
This part presents the calculation of the deterministic component of the direct solar radiation.

Calculation of the Geometrical Parameters
Before presenting the direct solar radiation by a sophisticated model, it is appropriate to define the geometrical parameters describing the Sun position in a point of the Earth at a given moment. To reach this objective, two coordinate systems are used; the first is the equatorial coordinate system which defines the position of the Sun's center relative to the Earth's center. It is presented by the Sun's declination (δ) and the hour angle ( ). The second coordinate system is the local coordinate system which defines the position of the Sun's center relative to the point of measurement. It is presented by the sun height (α) and the azimuth angle (ψ).
This section focuses on the definition of all the geometric parameters required to calculate the direct solar radiation on a horizontal surface [33][34][35][36][37][38][39][40][41][42][43][44]. It should be noted that the literature contains more than one expression for some parameters. Therefore, a comparative study has been performed, in this section, in order to choose the adequate equations for three parameters: the Sun's declination (δ), the equation of time (E) and the Earth-Sun distance correction factor K D . The Sun's declination (δ) is the first equatorial coordinate. It is the angle formed by the equator plane and the Earth-Sun direction ( Figure 4). The Sun's declination is approximated by Cooper as follows [33]: where, J is the rank of day in the current year (1 for 1 January). This expression calculates (δ) in degree, and the error is included in sin( ) = 0.397744 × sin( ) L is the true longitude, or ecliptic longitude of the Sun in degrees. Its expression is: C is the equation of the center, its expression is: Ma is the mean anomaly, and its expression is: = 357.5291 + 0.98560028 × (9) e = 0.1671, is the eccentricity of the ellipse.
To validate the Equations (5) and (6), their errors are drawn (Error 1 and Error 2 respectively) based on the IMCCE calculation in Figure.
The parameter 365.24 (in days) is an approximated value of the topic year. N is the rank of the day beginning on 1 January 2013.
The second expression involves more geometric parameters. It is defined as follows [33][34][35][36]: L is the true longitude, or ecliptic longitude of the Sun in degrees. Its expression is: C is the equation of the center, its expression is: Ma is the mean anomaly, and its expression is: e = 0.1671, is the eccentricity of the ellipse. To validate the Equations (5) and (6), their errors are drawn (Error 1 and Error 2 respectively) based on the IMCCE calculation in Figure 5. Error 1 and Error 2 don't exceed 0.3475 • . The mean value of Error 1 is 0.1514 • and of Error 2 is 0.2270 • for the year 2016. So the adopted equation is the Equation (5).   (5) and (6) and IMMCE calculation to choose the Sun's declination equation. Figure 5. Comparison between Equations (5) and (6) and IMMCE calculation to choose the Sun's declination equation.

Equation of Time
The period of Earth's sidereal rotation is around 23 h 56 min 4 s (=23.93 h), but the Sun returns to the Earth local meridian in 24 h on average. In fact, when the Earth turns around its axis, it moves on its orbit by approximately (360.98/365.25 = 0.98 • ). To make a rotation of 360.98 • it takes (360.98/360 × 23.93 = 24.00) h. It would be simple if the orbital speed was constant, and the axis of Earth rotation was orthogonal to the ecliptic plane, but the orbital speed of Earth is not regular. It is faster when the Earth approaches the Sun. Moreover the rotation axis is tilted relative to the normal direction of the trajectory plane. The angle that the Earth forms in relation with the Sun in the followed trajectory during the year varies from one year to another one, around one degree. Therefore the true time moves away periodically from the mean time: this gap is given by the equation of time E according to the day of the year. In a first order, E depends on the orbital eccentricity and the tilt of the rotation axis. This dependency is at the origin of the traditional expression which is the sum of three sinusoidal functions: where: The mean error of this equation is 25 s, and its maximal error is more than one minute. For more precision, the IMCCE astronomers have obtained an equation which takes into account the main global disturbances for the 1900-2100 period [39]. This equation contains 12 sinusoidal functions and 2 pseudo-periodic functions. But it can be curve fitted to obtain two other simplified equations for the 2013-2023 period [11]. The first equation is: The second expression is based on the equation of the center (C) and the influence of obliquity (R) [13,15]. It is: where, R is called the influence of obliquity. It is defined as: where, y = tan(ε/2) and ε is the tilt of the Earth axis, ε = 23.43.
To validate the Equations (12) and (13), their errors are drawn (Error 1 and Error 2 respectively) based on the IMCCE calculation in Figure 6. The maximal value of Error 1 is 14.77 s; its mean value is 8.16 s. However the maximal value of Error 2 is 5.598 s and its mean value is 2.532 s. As a result the adopted expression is the Expression (13).

Local and Solar Time
The local time is the civil time. It is linked to a reference time zone. Each territory is affiliated with one of the 24 time zone, depending on the geographic and politic considerations. The zero's time zone corresponds to the space situated between 7.5° E and 7.5° W from the meridian of origin, passing through London-Greenwich.
Within a time zone, such location can be situated between 7.5° E (+1/2 h) and 7.5° W (−1/2 h) from the central meridian of the time zone. This offset must be taken into account by the geometric correction term (Lon/15) (difference of 4 min in degree of longitude).
Lon is the longitude (Figure 7). Finally, the true solar time is deduced from the local time as following:

Local and Solar Time
The local time is the civil time. It is linked to a reference time zone. Each territory is affiliated with one of the 24 time zone, depending on the geographic and politic considerations. The zero's time zone corresponds to the space situated between 7.5 • E and 7.5 • W from the meridian of origin, passing through London-Greenwich.
Within a time zone, such location can be situated between 7.5 • E (+1/2 h) and 7.5 • W (−1/2 h) from the central meridian of the time zone. This offset must be taken into account by the geometric correction term (Lon/15) (difference of 4 min in degree of longitude).
Lon is the longitude (Figure 7). Finally, the true solar time is deduced from the local time as following: TCF is the civil time and cc is the time difference comparing to GMT. The sign of (Lon/15) is taken negative in the East of Greenwich and positive in its West.
Energies 2018, 11,620 zone corresponds to the space situated between 7.5° E and 7.5° W from the meridian of origin, passing through London-Greenwich. Within a time zone, such location can be situated between 7.5° E (+1/2 h) and 7.5° W (−1/2 h) from the central meridian of the time zone. This offset must be taken into account by the geometric correction term (Lon/15) (difference of 4 min in degree of longitude).
Lon is the longitude (Figure 7). Finally, the true solar time is deduced from the local time as following: TCF is the civil time and cc is the time difference comparing to GMT. The sign of (Lon/15) is taken negative in the East of Greenwich and positive in its West.

Hour Angle
The hour angle (ϖ) is the second equatorial coordinate. It is defined as the angle composed by the local meridian plane and the meridian plane of the Sun center. It is taken positive in the direction of the East (Figure 4). The hour angle expression is as follows:

Hour Angle
The hour angle ( ) is the second equatorial coordinate. It is defined as the angle composed by the local meridian plane and the meridian plane of the Sun center. It is taken positive in the direction of the East (Figure 4). The hour angle expression is as follows:

Sun Height
The Sun height (α), or its complementary, the solar zenith angle (θz), (θz = 90 • − α), is the first local coordinate. It is the angle between the line formed by the center of sun and the point of measurement, and the horizontal plane of the observer (Figure 8). The Sun height is defined as: sin(α) = sin(δ)· sin(ϕ) + cos(δ)· cos(ϕ)· cos( ) (17) (φ) is the latitude in degree (Figure 7). The solar zenith angle (θz) is the angle between the line formed by the center of Sun and the point of measurement, and the vertical plane of the observer (Figure 8).    The azimuth angle (ψ) is the second Sun local coordinate. It is the angle composed by the projection of the Sun direction in the horizontal plane and the South direction. The azimuth angle (ψ) is defined as follows: cos(ψ) = cos(δ)· cos( )· sin(ϕ) − sin(δ)· cos(ϕ) cos(α)

Solar Radiation at the Top of Atmosphere
The solar radiation at the top of the earth atmosphere (G 0 ) can be easily deduced from the previous equations. It is defined as follows: where E sc is the solar constant and its value is equal to 1367+/−0.1 W/m 2 [40]. It is the quantity of solar energy that receives a surface of 1 m 2 located at a distance of 1 au (astronomical unit) and exposed perpendicularly to the solar radiation, without considering the atmosphere. R m is the mean Earth-Sun distance (in astronomical units) and R(J) is the mean distance for the day J. The expression ( R m R(J) ) 2 is called the Earth-Sun distance correction factor and it is named K D .
Daniel and Gautret evaluate K D by the following equation (the error is around 0.2%) [41]: Because the Sun's declination and the Earth-Sun distance are linked (precession and nutation), K D can also be approached by an expression based on the Sun's declination and developed as follows [42]: The error of this equation is less than 0.15%. Besides, another more precise equation is proposed in [43]: Its error (error 1) is drawn in Figure 9. It does not exceed 0.11%. However the sum of a constant and a sinusoid cannot describe K D perfectly. This factor can be calculated thanks to the fact the Earth-sun Distance calculated by resolving the orbit equations. Based on the IMCCE calculation, a precise K D expression is obtained with a maximal error equal to 0.046% (Error 2 in Figure 9): As a result, the adopted equation is Equation (24). Using the chosen equations (sun's declination, time equation and distance correction factor earth-sun), G 0 (Equation (20)) is calculated and validated with National Aeronautics and Space Administration (NASA) calculations [44]. An error equal to 0.13 % is obtained for the year 2016, when the studied location is Biarritz (France).
(Error 2 in Figure 9): As a result, the adopted equation is Equation (24).  (23) and (24) and IMMCE calculation to choose the Earth-Sun distance correction factor.
Using the chosen equations (sun's declination, time equation and distance correction factor earth-sun), G0 (Equation (20)) is calculated and validated with National Aeronautics and Space Administration (NASA) calculations [44]. An error equal to 0.13 % is obtained for the year 2016, when the studied location is Biarritz (France).

The Direct Solar Radiation Model
Before modeling the "normal" sky it is interesting to model the sky without cloud cover. This concept is called the "clear sky" solar radiation model [45][46][47]. The estimation of the solar radiation by the clear sky model allows obtaining correct results when the sky is clear, particularly in summer [22,48,]. Several models give possible ways to represent the solar radiation through the clear sky,  (23) and (24) and IMMCE calculation to choose the Earth-Sun distance correction factor.

The Direct Solar Radiation Model
Before modeling the "normal" sky it is interesting to model the sky without cloud cover. This concept is called the "clear sky" solar radiation model [45][46][47]. The estimation of the solar radiation by the clear sky model allows obtaining correct results when the sky is clear, particularly in summer [22,48]. Several models give possible ways to represent the solar radiation through the clear sky, including Kasten model [49], Molineaux model [50], and Solis model [51]. In our study, the SOLIS model is used to describe the deterministic component of direct solar radiation. This model is a result of the European Héliosat-3 project of Oldenburg University. It gives excellent results when it is compared with measures obtained in Europe [22,48]. The simplified SOLIS model is an approximation of the Radiative Transfer Model (RTM) based on Lambert-Beer law [52].
The proposed direct solar radiation expression is as follows: ) (25) where, b is a constant adjustment parameter and τ is the optical depth.

Interpolation of Downloaded Data
Two interpolations are applied to the downloaded data in order to obtain the cloud cover in the exact location of ESTIA (latitude = 43.44, longitude = 1.55) with a time interval of 5 min (the time interval used for the weather station measurements). The first is 2-D geographical interpolation (related to longitude and latitude), and the second one is 1-D time interpolation. In order to provide best results, the piecewise polynomial interpolation [53,54] will be used. Three interpolation orders (from 1 to 3) are studied.

Geographical Interpolation
The first tests are performed using the first-order interpolation, named also linear interpolation (see Figure 10). The result is a set of quadrilateral zones spread on 0.25 • of longitude and latitude.
All the values are located between 0% and 100%. The inconvenience of the linear interpolation is that the obtained zones make ridges unfavorable for the training phase of the neural network. In order to avoid this disadvantage, the order of the interpolation is increased to two (spline) and three (cubic). As it can be seen in Figure 11, the surface obtained by the interpolation is smooth. The inconvenient of spline and cubic interpolations is that some extrema are not located between 0% and 100%. To remedy Energies 2018, 11, 620 13 of 21 that situation, the solution is to limit the values to the interval [0,100]. The exceeding of the limit interval is smaller for the cubic interpolation, so it will be used thereafter.
The first tests are performed using the first-order interpolation, named also linear interpolation (see Figure 10). The result is a set of quadrilateral zones spread on 0.25° of longitude and latitude. All the values are located between 0% and 100%. The inconvenience of the linear interpolation is that the obtained zones make ridges unfavorable for the training phase of the neural network. In order to avoid this disadvantage, the order of the interpolation is increased to two (spline) and three (cubic). As it can be seen in Figure 11, the surface obtained by the interpolation is smooth. The inconvenient of spline and cubic interpolations is that some extrema are not located between 0 and 100%. To remedy that situation, the solution is to limit the values to the interval [0,100]. The exceeding of the limit interval is smaller for the cubic interpolation, so it will be used thereafter.

Time Interpolation
For each time step (3 h) of the geographic interpolation we obtain a value corresponding to the sailboat position. The time interpolation is used to compute values for each 5 min. When the time interpolation is linear, the result is a set of line segments spread on 3 h intervals. In this case, the transition between line segments is sudden and this is not favorable for the neural network training. Using a spline interpolation, the transition between the downloaded data is harmonious. According to the data exceeding of the interval [0,100], a limitation of the output must be applied. It should be noted that this exceeding don't occur frequently in the case of time interpolation, so the order is maintained to two (spline).

Time Interpolation
For each time step (3 h) of the geographic interpolation we obtain a value corresponding to the sailboat position. The time interpolation is used to compute values for each 5 min. When the time interpolation is linear, the result is a set of line segments spread on 3 h intervals. In this case, the transition between line segments is sudden and this is not favorable for the neural network training. Using a spline interpolation, the transition between the downloaded data is harmonious. According to the data exceeding of the interval [0,100], a limitation of the output must be applied. It should be noted that this exceeding don't occur frequently in the case of time interpolation, so the order is maintained to two (spline).

Evaluation Criteria
The goal of this research study is to estimate the solar energy availability. On the assumption that the energy is determined as the integral of the power, we can consider that the most important validation criteria is the Daily Mean of the Power Error DMPE (in W/m 2 ). The DMPE identifies the daily excess or the daily lack of the prediction results compared to the real measurements. The benefit of the DMPE use is that its calculation is performed to ensure that the errors having an opposite signs cancel each other, which matches with the objective of this research. The DMPE is defined as follows: where N denotes the number of pattern pairs, y i andŷ i are the measured and predicted solar radiation of the i th pattern pair respectively.
Since the DMPE is difficult to integrate into the training process, the Mean Square Error (MSE) is used in this phase and maintained in the phase of simulation results evaluation. The MSE is defined as follows: It should be noted that, the MSE is calculated for normalized data (0.05 to 0.95) and the DMPE is calculated for actual size data (0 to around 1000 W/m 2 ).

Results and Discussion
Several tests are considered to select the database structure and the neural network configuration. Before starting simulations, the dataset is separated into three parts: training, testing and validation sets. The training phase uses the training set to compute the weights and bias of the neural network and the test set to test it. After finishing this phase the validation set is used to simulate the model and evaluate its performance. All simulations are performed using normalized dataset in order to adjust data and not saturate the neurons.

Choice of the Dataset Structure
The first simulations of this research study are performed using a set of test and training containing several days of the year. Actually, the idea is to use a database which can represent all the year in order to obtain a sufficient generalization capability and facilitate the training process. It is assumed that the use of a dataset containing the 365 days of the year is not beneficial for several reasons. The use of more than a certain threshold of data increases the calculation time and mays cause an overtraining. In addition some extreme phenomena can never be fully predicted even if the dataset is extended. As a result, a chosen number of days is selected to perform the training and test phase, as analyzed later.
As it can be seen in Figure 12, the neural network cannot follow correctly the direct solar radiation curve, particularly at night when the solar radiation must be null. The observed leaps in the predicted curve are explained by the fact that the set of test and training consists of data of non-successive days. Indeed, the use of these data generates two forms of inconsistency: the cloud cover inconsistency and the solar characteristics inconsistency. The cloud cover inconsistency is observed in the moment of the transition between non-successive days (at midnight) because the meteorological parameters are different. This transition engenders a linearity problem for the neural network training. The solar characteristics inconsistency is due to the temporal gap between data, which causes changing features between days such as the daily solar peak and the daily sunshine hours. This gap engenders a generalization problem for the neural network. On the premise of the daily training, a compromise must be found in order to obtain the adequate database size which allows to obtain good results without saturating the calculation process. Several training are performed and the best results are briefly presented in Table 2. A test and training dataset of 10 days allows to obtain the minimal error. In this case, the MSE and the DMPE are 0.00695 and 41.19645 W/m 2 respectively. The global solar radiation of the weather station is measured with a time step of 5 min. Therefore the interpolated cloud cover and the calculated direct solar radiation have the same time step. However, as observed in the training with the previous dataset (beginning of 6.2.1), the output of the neural network is not meant to follow the rapid fluctuations of the measured solar radiation ( Figure  12). Applying a filter on the input signal, it becomes more pertinent because its characteristic is close to the output signal. Thus it is interesting to consider the average of the measured data in the prediction process. Therefore, several averages are performed and the best prediction results are stored in Table 3. There are moving averages of 10 and 30 min over intervals of one hour. It may be concluded that the two simulation results are similar. Therefore, in order to minimize errors and computing time, the chosen interval average is 30 min. The second important step in the determination of the solar radiation predictor is to find the adequate neural network structure. Several simulations are performed in order to select the different parameters of the NARX neural network model. This section cannot contain all detailed results. Therefore some parameters choices are presented in Table 4, and the decisions related to the three most important parameters are described: the choice of the activation function in each layer, the choice of the neuron number in each layer and the use of weights in the daily trainings. For these reasons, the idea is to use a set of test and training data containing a number of consecutive days and to repeat the training process one time per day to forecast solar radiation for one coming day. Using this approach, the used database becomes pertinent because it consists of a number of days just before the day concerned by the prediction. However, a daily training consumes calculation resources and can engender discontinuities if the radiation is not equal to zero at the moment of the training. Therefore, the choice of the training time is important. Since the solar radiation is null at night, and the training should be achieved every day, we have chosen to perform it at midnight, the transition moment between days.
It must be noted that the daily training is useful for the solar radiation prediction in the sailboat. In fact, the sailboat provokes the variation of two important parameters:

•
The meteorological conditions influencing the cloud cover variations.

•
The localization influencing the direct solar radiation model variations, especially as in this study the localization data (longitude and latitude) are intercalated in a preliminary calculation (the SOLIS model) and not as inputs of the neural network.
On the premise of the daily training, a compromise must be found in order to obtain the adequate database size which allows to obtain good results without saturating the calculation process. Several training are performed and the best results are briefly presented in Table 2. A test and training dataset of 10 days allows to obtain the minimal error. In this case, the MSE and the DMPE are 0.00695 and 41.19645 W/m 2 respectively. The global solar radiation of the weather station is measured with a time step of 5 min. Therefore the interpolated cloud cover and the calculated direct solar radiation have the same time step. However, as observed in the training with the previous dataset (beginning of Section 6.2.1), the output of the neural network is not meant to follow the rapid fluctuations of the measured solar radiation (Figure 12). Applying a filter on the input signal, it becomes more pertinent because its characteristic is close to the output signal. Thus it is interesting to consider the average of the measured data in the prediction process. Therefore, several averages are performed and the best prediction results are stored in Table 3. There are moving averages of 10 and 30 min over intervals of one hour. It may be concluded that the two simulation results are similar. Therefore, in order to minimize errors and computing time, the chosen interval average is 30 min. The second important step in the determination of the solar radiation predictor is to find the adequate neural network structure. Several simulations are performed in order to select the different parameters of the NARX neural network model. This section cannot contain all detailed results. Therefore some parameters choices are presented in Table 4, and the decisions related to the three most important parameters are described: the choice of the activation function in each layer, the choice of the neuron number in each layer and the use of weights in the daily trainings. The choice of the activation function in each layer of the neural network proves to be an important constituent element. The activation function used in the neurons of the input and hidden layer is the sigmoid. This function is especially advantageous in the neural networks which are trained by back-propagation algorithms (Levenberg-Marquardt algorithm). Furthermore, the sigmoid function presents the advantage to be derivable, which makes the weights learning of the neural network easier.
The activation functions of the output layer in the majority of neural studies are linear. However, this choice may change depending on the application. Figure 13 shows the training result using linear function in the output layer. The trained curve is able to follow the real curve from sunsets to sunrises, but it is not able to simulate the solar radiation behavior during daylight hours. When the sigmoid or tansig is used, the neural network is able to approximate the real solar radiation curve in the daytime. The sigmoid is unable to stabilize the output signal to zero at night, because the function codomain is ]0, 1[. The codomain of the tansig function is ]−1, 1[ and the behavior of the tansig is similar to the linear function around zero, thus this value is easily reachable by the output signal. Furthermore, the tansig function in the output layer has proved its performance in the generalized MLP architectures [55]. For these reasons, it will be used thereafter.
The sigmoid is unable to stabilize the output signal to zero at night, because the function codomain is ]0, 1[. The codomain of the tansig function is ]−1, 1[ and the behavior of the tansig is similar to the linear function around zero, thus this value is easily reachable by the output signal. Furthermore, the tansig function in the output layer has proved its performance in the generalized MLP architectures [55]. For these reasons, it will be used thereafter. In order to complete the structure of the NARX neural network model, a last important parameter must be studied: the neuron number in each layer. Table 5 presents the best results obtained when the number of neurons is varied in each layer of the NARX neural network.  In order to complete the structure of the NARX neural network model, a last important parameter must be studied: the neuron number in each layer. Table 5 presents the best results obtained when the number of neurons is varied in each layer of the NARX neural network.  Table 5 shows that the best structure of the neural network is obtained using 15 neurons in the input layer, 15 neurons in the hidden layer and one neuron in the output layer. The MSE and the DMPE obtained with this structure are 0.00410 and 30.4164 W/m 2 respectively.
Another issue concerns the initialization of the ANN weights and biases. In fact, normally they are randomly generated. But on the assumption that the training is periodic, one idea can be to save the weights and biases vectors of the first training (which predicts the solar radiation for the first day) and to use it for the coming days. Therefore two different ways of ANN calibration are tested and compared:

•
The weights and biases vectors are randomly generated only once, in the first training phase. Then, in the following periodic training phases, the same weights and biases values are used. • Weights and biases are randomly initialized by ANN in each periodic training phase.
As shown in Table 6, best performances are achieved when ANN generate weights and biases vectors randomly for each training phase. In conclusion, the best obtained error performance is 0.00279 for MSE and 24.0584 W/m 2 for DMPE. Figure 14 shows an example of predicted direct solar radiation day using the obtained NARX neural network model. The predicted curve doesn't follow the rapid fluctuations of the solar radiation real curve but it follows the general curve and it leads to good results, in particular when considering the DMPE. In conclusion, the best obtained error performance is 0.00279 for MSE and 24.0584 W/m 2 for DMPE. Figure 14 shows an example of predicted direct solar radiation day using the obtained NARX neural network model. The predicted curve doesn't follow the rapid fluctuations of the solar radiation real curve but it follows the general curve and it leads to good results, in particular when considering the DMPE.

Conclusions
This paper proposes a NARX neural network model for the direct solar radiation prediction on a horizontal surface. The main finding of this research study is that the training phase of the neural Figure 14. Predicted curve of the direct solar radiation when using the best obtained simulation results.

Conclusions
This paper proposes a NARX neural network model for the direct solar radiation prediction on a horizontal surface. The main finding of this research study is that the training phase of the neural network is performed periodically, taking into consideration several parameters, such as the cloud cover and solar characteristics and the sailboat mobility.
Several simulations carried out with different evaluation criteria are performed and evaluated using MSE and DMPE errors. The best results (0.00279 for MSE and 24.0584 W/m 2 for DMPE) are obtained for: a dataset of 10 days with a moving average of 30 min over intervals of one hour, a NARX model consisting of 15 neurons on the input and the hidden layers, an input and hidden layer which contain the sigmoid function, an output layer which contain the tansig function, and a random initialization of weights.
The developed predictor is useful to obtain the direct solar radiation on a mobile horizontal surface, but it does not consider the two tilt angles needed to take into account the sailboat pitch and roll movements. In addition, the shadow of the sails is not considered either. Therefore, the designed predictor will be improved to consider these new constraints. The new predictor will be developed in the frame of the aforementioned project about a 100% renewable race sailboat, first to predict the direct solar radiation on a sloped and potentially shady surface, and then the amount of available power from PV arrays installed in the boat. This new predictor will be used by the energy management system of the boat.