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

Abstract

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, artificial 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 specific 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 specific 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.

1. 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 Autoregressive Exogenous (NARX) model. In [6], the authors use only the historical solar radiation as input of the forecaster, so the wavelet transformation is used to split the historical data into adapted series for the forecast. The authors apply the ARMA concept as a linear forecaster and the NARX model as a nonlinear form recognition instrument, in order to decrease the errors of the wavelet-ARMA forecast.

In [2] also, only the historical data of the solar radiation is used and decomposed. This decomposition is performed via Local Mean Decomposition (LMD) and Empirical Mode Decomposition (EMD). The result is a set of nonlinear sub-components, so the authors of [2] combine two nonlinear prediction models. They use Least Squares Support Vector Machine (LSSVM) and Volterra models to predict the high-frequency subcomponents and low-frequency subcomponents, respectively. The final result is obtained by superimposing the sub-prediction results.

In [7] four models of machine learning are compared: Random Forests (RF), neural network based on Multi-Layer Perceptron (MLP), Linear Regression (LR) and k-Nearest Neighbors (kNN). The objective is to demonstrate the efficiency of nonlinear models to forecast two times series: the hourly solar radiation on horizontal surface and the solar irradiance obtained on PV panels at different tilted surfaces. The forecast of the solar radiation is performed using three types of data inputs: ground weather data, satellite remote sensing data and calculations of sun position parameters (declination, hour, zenith, elevation and azimuth angles).

References [3,4,5,8] forecast the daily solar radiation using ANNs. In [3], its authors use the correlation criteria to determine the endogenous and the exogenous inputs to take into account. Four endogenous time lags are taken for the clear sky model of solar radiation. From several meteorological parameters only three are selected as exogenous inputs: relative humidity, sunshine duration and cloud cover.

Based on the nature of the obtained data, the authors of [4] divide data differently: predicted data and statistical data, so they combine: 1/several statistical data as the calculated cloud ratio, the maximum hourly variation of the solar radiation, the solar radiation absolute daily variation between the day (t) and the day (t − 1), etc., and 2/numerical weather prediction data which contain the one day ahead forecast of the cumulative radiation each 3 h. As a result, a hybrid dataset is used as neural network input.

In [5], the authors don’t classify the input data, but rather they train four neural networks with four combinations of input features, with the aim of considering the influence of several meteorological parameters on the prediction results. They conclude that the best performance is obtained when using the following inputs: the year day, the mean daily solar radiation at the top of atmosphere, the maximum sunshine hours, the mean daily air temperature, the mean daily relative humidity and the wind speed.

In some other studies, the neural network technique is combined with other nonlinear models. For example, in [8], the authors propose a diagonal recurrent wavelet neural network, in order to estimate the hourly solar radiation of the next day. It is a model founded on the combination of wavelet neural network with recurrent ANN and fuzzy technology, where the wavelet basis is implemented as activation function for the neurons and the input vector of the neural network contains defuzzificated data of the nebulosity.

The results of these researches are advanced and give much better performances than the conventional prediction models. However, they cannot be generalized. They don’t satisfy the requirements of certain specific situations. Among these situations, that of this research, whose main goal is the prediction of the direct solar radiation on a horizontal surface, in order to predict the solar radiation on a sloped surface, then the quantity of obtainable power from PV arrays in a race sailboat. Contrary to the aforementioned research studies, the proposed solution in this paper takes the change of localization into consideration and allows to take into account the other specific conditions (tilted surface, shadow …) of the sailboat operation that will be considered in our future studies.

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.

2. 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 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:

$$\widehat{y}\left(t+1\right)=F\left(\begin{array}{l}y\left(t\right),y\left(t-1\right),\dots ,y\left(t-n_y\right),x\left(t+1\right),\\ x\left(t\right),x\left(t-1\right),\dots ,x\left(t-n_x\right)\end{array}\right)$$

(1)

$$\widehat{y}\left(t+1\right)=F\left(\begin{array}{l}\widehat{y}\left(t\right),\widehat{y}\left(t-1\right),\dots ,\widehat{y}\left(t-n_y\right),x\left(t+1\right),\\ x\left(t\right),x\left(t-1\right),\dots ,x\left(t-n_x\right)\end{array}\right)$$

(2)

where, F(·) is the mapping function of the neural network, $\widehat{y}$(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). $\widehat{y}$(t), $\widehat{y}$(t − 1), …, $\widehat{y}$(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.

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 $\widehat{y}$(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\times w_{ij}$. An activation function f, given by

Table 1. Examples of activation functions.

Function	Definition	Codomain
Linear	x	The Real domain IR
Sigmoid (Logsig)	$\frac{1}{1+e^{-x}}$	]0, 1[
Hyperbolic tangent (Tansig)	$\frac{e^x-e^{-x}}{e^x+e^{-x}}$	]-1, 1[

, is then performed to obtain the following neuron output:

$$y_i=f({{\displaystyle \sum }}_{j=1}^n{x}_j·w_{ij})$$

(3)

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.

3. 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 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.

4. The Deterministic Component of the Direct Solar Radiation

This part presents the calculation of the deterministic component of the direct solar radiation.

4.1. 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.

4.1.1. 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]:

$$\delta =23.45·\sin (\frac{2\pi (J+284)}{365})$$

(4)

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 [−1.4°, +0.5°].

For more precision, two other expressions can be obtained by the curve fitting of data calculated by the “Bureau des Longitudes” for the period 2013–2023 (Institut de Mécanique Céleste et de Calcul des Éphémérides (IMCCE)) [34]. The first expression is:

$$\delta =0.38+23.26·\sin (\frac{2\pi N}{365.24}-1.395)+0.375·\sin (\frac{4\pi N}{365.24}-1.47)$$

(5)

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]:

$$\sin \left(\delta \right)=0.397744\times \sin \left(L\right)$$

(6)

L is the true longitude, or ecliptic longitude of the Sun in degrees. Its expression is:

$$L=280.4665+C+0.98564736\times N$$

(7)

C is the equation of the center, its expression is:

$$C=\frac{180}{\pi }[\left(2e-\frac{1}{4}{e}^3\right)\sin \left(Ma\right)+\frac{5}{4}{e}^2·\sin \left(2Ma\right)+\frac{13}{12}{e}^3·\sin \left(3Ma\right)]$$

(8)

Ma is the mean anomaly, and its expression is:

$$Ma=357.5291+0.98560028\times N$$

(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 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).

4.1.2. 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:

$$E=7.53·\cos \left(B\right)+1.5·\sin \left(B\right)-9.87·\sin \left(2B\right)$$

(10)

where:

$$B=\frac{2\pi \left(J-81\right)}{365}$$

(11)

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:

$$E=7.36·\sin \left(\frac{2\pi N}{365.242}-0.071\right)+9.92·\sin \left(\frac{4\pi N}{365.242}+0.357\right)+0.305·\sin \left(\frac{6\pi N}{365.242}+0.256\right)$$

(12)

The second expression is based on the equation of the center (C) and the influence of obliquity (R) [13,15]. It is:

$$E=\left(C+R\right)·4$$

(13)

where, R is called the influence of obliquity. It is defined as:

$$R=\frac{180}{\pi }·(-y^2·\sin (2L)+\frac{y^4}{2}·\sin (4L)-\frac{y^6}{3}·\sin (6L))$$

(14)

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).

4.1.3. 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:

$$TS=TCF-cc+\frac{E}{60}\pm \frac{\mathrm{Lon}}{15}$$

(15)

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.

4.1.4. 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:

$$\varpi =15·\left(TS-12\right)$$

(16)

4.1.5. 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 \left(\alpha \right)=\sin \left(\delta \right)·\sin \left(\phi \right)+\cos \left(\delta \right)·\cos \left(\phi \right)·\cos \left(\varpi \right)$$

(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).

4.1.6. Azimuth Angle

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:

$$\sin \left(\psi \right)=\frac{\sin \left(\delta \right)·\sin \left(\varpi \right)}{\cos \left(\alpha \right)}$$

(18)

$$\cos \left(\psi \right)=\frac{\cos \left(\delta \right)·\cos \left(\varpi \right)·\sin \left(\phi \right)-\sin \left(\delta \right)·\cos \left(\phi \right)}{\cos \left(\alpha \right)}$$

(19)

4.1.7. Solar Radiation at the Top of Atmosphere

The solar radiation at the top of the earth atmosphere (G₀) can be easily deduced from the previous equations. It is defined as follows:

$$G_0={(\frac{R_m}{R(J)})}^2·E_{sc}·\sin (\alpha )$$

(20)

where E_sc is the solar constant and its value is equal to 1367+/−0.1 W/m² [40]. It is the quantity of solar energy that receives a surface of 1 m² 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 ${(\frac{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]:

$$K_D=1-0.034·\cos (\frac{2\pi J}{365})$$

(21)

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]:

$$K_D=1-0.034·\sin (\frac{2\pi (J-94)}{365})$$

(22)

The error of this equation is less than 0.15%.

Besides, another more precise equation is proposed in [43]:

$$K_D=1+0.0334 \cos \left(\frac{2\pi ·\left(N-2.7206\right)}{365.25}\right)$$

(23)

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):

$$\begin{array}{ll}{K}_D=1.000138& +0.03341·\cos \left[\frac{2\pi N}{365.2422}-0.051\right]+0.000699·\sin \left[\frac{4\pi N}{365.2422}+1.474\right]\\ & +0.000062 ·\sin \left[\frac{12.37\times 2\pi N}{365.2422}+2.2\right]\end{array}$$

(24)

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₀ (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).

4.2. 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:

$$G=G_0·\exp (-\frac{\tau }{{\sin }^b\left(\alpha \right)})$$

(25)

where, b is a constant adjustment parameter and τ is the optical depth.

5. 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.

5.1. 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 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.

5.2. 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).

6. Adjustments, Results and Discussion

6.1. 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²). 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:

$$\mathrm{DMPE}=\frac{{\displaystyle {\sum }_1^N\left(y_i-{\widehat{y}}_i\right)}}{N}$$

(26)

where N denotes the number of pattern pairs, y_i and $\widehat{y}$_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:

$$\mathrm{MSE}=\frac{{\displaystyle {\sum }_1^N{\left(y_i-{\widehat{y}}_i\right)}^2}}{N}$$

(27)

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²).

6.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.

6.2.1. 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.

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. Choice of the test and training size of the database.

Database Size	MSE	DMPE (W/m2)
5 days	0.011	60.228825
10 days	0.00695	41.19645
15 days	0.009625	50.2089

. 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² 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. Choice of the time average of the database.

Time Average	MSE	DMPE (W/m2)
10 min	0.00732	44.4344
30 min	0.00695	41.19645

. 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.

6.2.2. Choice of the Neural Network Structure

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. Parameters of the NARX model structure.

Property	Choice
Number of hidden layers	1
Normalization Interval of dataset	[0.05, 0.95]
Delay vectors	Input data: [0, 1] Target: [1, 2]
Training parameters	Error: MSE Learning algorithm: Levenberg-Marquardt

, 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.

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. Choice of the neurons number in the different layers.

Number of Neurons	MSE	DMPE (W/m2)
10 × 10 × 1	0.00724	59.5724
15 × 15 × 1	0.00410	30.4164
16 × 16 × 1	0.01438	73.4646
20 × 20 × 1	0.00768	45.0513
22 × 22 × 1	0.00695	41.1964

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² 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. Choice of registration of weights for periodic trainings.

Generation of Weights	MSE	DMPE (W/m2)
Initialized	0.00410	30.4164
Random	0.00279	24.0584

, 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² 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.

7. 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² 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.