A Forecasting-Based Control Algorithm for Improving Energy Managment in High Concentrator Photovoltaic Power Plant Integrated with Energy Storage Systems

The High Concentrator Photovoltaic (HCPV) technology, due to its high efficiency, is considered one of the most promising solutions for the exploitation of sun-irradiation-based Renewable Energy Sources (RES). Nevertheless, the HCPV production is strictly connected to the Direct Normal Irradiation (DNI) making this photovoltaic technology more sensible to cloudiness than traditional ones. In order to mitigate the power intermittence and improve production programmability, the integration between Energy Storage Systems (ESSs) and HCPV, resorting to forecasting algorithms, has been investigated. Specifically, a local weather forecasting algorithm has been used for estimating the daily time evolution of DNI, air Temperature (T), Wind Speed (WS), and Air Mass (AM). These data are subsequently processed by means of an accurate HCPV model for the estimation of one day-ahead daily power production profile. The processing of HCPV forecasted generation by means of a properly tuned filter-based algorithm allows one day-ahead the definition of power profiles of ESS and power plant respectively, considering also the ESS constraints and the characteristic of the implemented real-time control algorithm. The effectiveness of the proposed forecasting model and control algorithm is verified through a simulation study referring to the solar power plant constituted by HCPV and ESS installed in Ottana, Italy. The results highlight that the application of the proposed approach lessens the power fluctuation effect caused by HCPV generation preserving the batteries at the same time. The feasibility and advantages of the proposed approach are finally presented.


Introduction
Nowadays, photovoltaic can be considered a mature technology and it will certainly have in the next future a fundamental role in electricity generation from Renewable Energy Sources (RES). In particular, the High Concentrator Photovoltaic System (HCPV) stands out among the available photovoltaic technologies thanks to the high efficiencies it is characterized by. Compared to the silicon-based monocrystalline cells, the average efficiency of HCPV cell is two times higher, reaching values equal to 41% for commercial cells [1], up to 44% in laboratories [2,3]. Such high values are obtainable thanks to the multi-junction technologies that consist of using multiple layers of different materials each of which reacts to different sunlight wavelengths. This kind of cell achieves its best performance for a specific level of irradiation. For this purpose, an optical concentrator that collects

The Case of Study
The Ottana solar facility is sited in Sardinia, Italy. Its main aim is to evaluate the reliability, potentiality, and effectiveness of Concentrating Solar Power (CSP) technologies for distributed generation applications characterized by power sizes lower than 5 MW. The Ottana concentrating solar power plant integrates a small scale CSP and a concentrating photovoltaic (CPV) system with electrochemical ESS and Thermal Energy Storage systems (TES) in order to evaluate the performance obtainable by concentrating solar facility in planning one-day ahead power output. In order to be able of managing the power fluctuations, the energy storage systems have been properly designed and controlled considering the support of weather forecasting service.
In particular, the proposed management strategy consists of two control algorithms: the first is addressed to synthesize the "one-day ahead" optimal scheduling, and the second performs a real-time smoothing control. The Ottana solar power plant is located in the center of Sardinia (Italy), and it is characterized by a 630 kWe CSP plant and a 400 kWe HCPV plant. As previously mentioned, the facility has been provided with both a thermal and an electric storage system. The CSP side is made up of a solar field based on Fresnel technology coupled to a power block based on Organic Rankine Cycle (ORC) unit. The ORC technology allows for operating with partial and with different input conditions without decreasing drastically its efficiency. Moreover, ORCs below 5 MWe are characterized by an interesting conversion efficiency, from thermal energy to electric energy.
The CSP field covers a net aperture area of about 8600 m 2 , and it generates a thermal power of about 5 MW at a reference irradiation condition. A diathermic oil works as heat transfer fluid and is characterized by a working temperature of 165 • C (inlet) and 275 • C (outlet). CSP is also equipped with two tanks for thermal storage, characterized by a capacity of about 15 MWhth, which means about 5 h of full-load operation. The HCPV side is constituted by 36 two-axis solar trackers, characterized by an overall power of about 400 kWp. A four-quadrant solar sensor and an onboard meteorological station, for measuring the DNI and the wind speed, have been installed on each solar tracker. The system combines astronomical ephemerides and solar tracker information in order to follow continuously the sun position.
The ESS consists of a Sodium-Nickel (NaNiCl) battery characterized by a capacity of 430 kWh and a maximum power of 300 kW. The ESS has been sized considering the expected power fluctuations occurring on HCPV systems as well as their dynamics, which results in being much greater if compared to those of the CSP plant. In particular, the system has been designed in order to allow the compensation of forecasting errors and a predictable generation of HCPV characterized by a time interval of 30 s.
The block diagram of the solar facility and the overview are reported in Figures 2 and 3, respectively.

HCPV Production Forecasting
The present section regards the forecast procedure. The first subsection presents the weather forecast algorithm, while the second one describes the five linear models adopted to forecast the power production on the basis of the forecast weather parameters.

Forecasting Method
Solar radiation forecasting can be achieved using different approaches, such as NWP, CMV, and time series analysis, the application and effectiveness of these techniques depend both on the weather-climate characteristics of the site and the lead time of the forecast itself.
DNI forecasting, in particular, is a challenge, given the greater variability with respect to e.g., wind intensity, temperature, or Global Horizonthal Irradiation (GHI), mainly due to fluctuations in cloud cover [11]. The research work is also affected by the low availability of DNI measurements.
In order to provide an operational service for DNI forecasting, for a specific site, it has been chosen in this case to exploit the potential of a non-specific weather forecast service (DarkSky https: //darksky.net). The service for the area of interest and the forecast horizon is mainly based on the results of National Oceanic and Atmospheric Administration (NOAA)'s Global Forecast System (GFS) global atmospheric circulation model (https://darksky.net/attribution) does not provide a direct estimate of GHI or DNI, but provides, in addition to air temperature and ground wind speed, an estimate of cloud cover, on an hourly basis and with a forecast horizon of one week. The cloud cover data are used to obtain an estimate of both GHI and DNI irradiance through a modulation of irradiance under Clear Sky conditions.
In particular, by following the approach described in [12], the GHI radiation can be calculated as: where GH I cs is the GHI in Clear Sky conditions calculated according to [13], and CC is the expected Cloud Cover, null in Clear Sky conditions and unitary for overcast sky.
The o f f set parameter allows for calibrating the relationship between CC and the ratio between GH I and GH I cs , in [12], which is taken equal to 0.35. DNI radiation is calculated in a similar way, using the approach described in [14] : where DN I cs represents the estimate of the Clear Sky direct radiation assessed according to [13].
The o f f set parameter in [14] is null. The Diffuse Horizonthal Irradiation (DHI) finally results: where z is the zenith angle of the Sun. Equation (1) was based on cloud cover forecasts obtained from two publicly-available NWP models: the North American Mesoscale Forecast System (NAM) of NOAA and the Regional Deterministic Prediction System (RDPS) of the Canadian Meteorological Center, while global radiation measurements were collected through two pyranometers located in California. Equation (2) was derived from the RDPS cloud cover forecast and direct irradiance data for stations located in the United States. It remains to be verified whether these models can be valid using forecasts from a different weather service and whether they are applicable to a different geographical area.

HCPV Production Linear Model
The models described in the following subsections belong to the family of models that evaluate HCPV production by means of a linear regression based on specific atmospheric parameters.
The models analyzed are those that best suit their application on the experimental power plant in Ottana. In fact, none of these require invasive procedures as the measurement of the cell temperature or the determination of the I-V characteristic, nor the application of algorithms based on neural networks. Instead, only the measurement of environmental and production parameters is needed, which is already implemented in the power plant and therefore available at no additional costs.
Three of the following models consider the AM as parameter. For its determination, the following equation by Sandia's laboratory has been used [15]: where z represents the zenith angle, p the local pressure in bar, and p 0 the reference pressure in bar. The value of the zenith angle (z) is determined by the procedure explained in [16]. For each model, the coefficient associated with a specific parameter is expressed as c n .

Model of the Standard ASTM
The ASTM model considers as parameters the DNI, the T, and the WS. The estimated ideal power production is expressed as: It is worth noting that the influence of the spectral effects is not taken into account and this is significant considering multijunction solar cell modules [8].

Model of E.F. FernáNdez et al.
The model developed by Fernández et al. [9] considers the spectral influence on the power production by using the AM. Thanks to this approximation, there is no need for any spectroradiometer or isotype cell [8]. The following model considers the influence of AM negligible when AM thre value is less than 2 [9]. Therefore, the characteristic equations are:

Modified Model of the Standard ASTM
In order to consider the AM influence, as explained above, in the equation of the ASTM, a term associated with the AM has been added. Therefore, the AM is characterized by the equation:

Modified Model of E.F. FernáNdez et al.
As Fernández explained in [9], the wind influence is negligible for wind speed values lower than 2 m/s. If the wind speed is higher than this value for almost the whole year, a term that considers its effect can be added to the equation of the Fernández model. In this model, AM thre is still considered equal to 2. Therefore, the equations related to the Modified-Fernández are as follows:

A Novel Hybrid Model
By the hybridization of the Modified ASTM and the Modified Fernández, a novel model is developed and presented in this paper. The proposed model considers all the previous shown parameters. In particular, both WS and AM are taken into consideration because it is not negligible, in particular considering windy sites (e.g., for the Ottana power plant location). This model is based on the Modified ASTM's equation but involves the effect of the AM, only if its value is less than AM thre , as in the Fernández model and in the Modified Fernández one. The equations associated with the Hybrid model are as follows: In this model, AM thre is set equal to 1.45.

Battery Energy Management
FBM is an energy management system strategy characterized by the decoupling between high frequency power profile evolution and the low frequency one, as described in [10]. The decoupling action is provided by a properly tuned low-pass filter. In particular, an FIR filter has been chosen. Since the filter acts on forecasting data (P hcpv ) and not on real-time samples, it has been updated as a zero-phase filter by applying it a second time on the flipped forecasting data array, allowing for avoiding any phase shift. In Figure 4a, an example of filter application is shown. The FIR output signal (P HCPV, f il ) is subsequently processed by the FBM algorithm. The ESS power setpoint returned by the FBM algorithm must observe the following power and energy constraints: P ess,min ≤ P ess (t) ≤ P ess,max e ess,min ≤ e ess (t) ≤ e ess,max (10) with P ess,min and P ess,max , maximum and minimum battery power values respectively; e ess,min and e ess,max are the minimum and maximum allowable battery energy values, respectively. Since a two threshold FBM algorithm has been chosen [10], the equation for calculating battery power (P ess ) is: P ess = −P HCPV, f il +r + 0 P ess,min + −P HCPV, f il +r − P ess,max 0 (11) in whichr + andr − are the upper and lower threshold, respectively; y x means "saturated between x and y". In this specific HCPV application, this equation has been modified as follows in order to compensate forecast power drops: P ess = −P HCPV, f il +r + 0 P ess,min + −P HCPV, f il +r − P ess,max 0 + P HCPV, f il − P hcpv P ess,max P ess,min The equation for battery energy calculation is as follows: with η d,ess and η d,ess as battery discharging and charging efficiency, respectively. Both thresholds are initialized to zero. For each iteration, the FBM algorithm calculates power and energy profiles; if one of the constraints in (10) is violated, it opportunely increases the upper threshold or decreases the lower one until no violations are observed. An additional modification of the algorithm proposed in [10] consists of a second application of the same algorithm that is executed subsequently in order to force the battery SoC, at the end of the day, to a chosen range (between 50% and 60% of SoC). This "refine" algorithm has the same logic of the main one, but it analyzes only the last part of the forecast array, from when e b reaches 60% of SoC to the end. In Figure 4b, an example of FBM application is shown, in which it is worth noting the peak shaving effect.

Smoothing Algorithms
HCPV energy production is characterized by high frequency fluctuation that can be even faster than traditional photovoltaic (PV) because of the strict dependence to DNI. The transmission system operators (TSO) have a limited response capacity and the greater the percentage of energy production from HCPV/PV, the more these fluctuations can lead to damaging system failures, above all in isolated cases such as islands. Coupling the HCPV/PV plant with an energy storage system with a proper energy management algorithm can lessen this kind of fluctuation. The smoothing algorithms meet these requirements by acting on power ramps inclination. Given any power profile, it is possible to compute the power variation associated with two consecutive values by the following equation [5]: with P Nom representing the nominal power of the plant. On the basis of this parameter, different smoothing algorithm can be developed. A generic block scheme of a smoothing algorithm is presented in Figure 5. It is worth noting that P ESS is defined as: where P g is the output of the smoothing algorithm and represents the power flowing to the grid. The non-ideal efficiency of the ESS has been considered introducing η pc and η bat that stand for power converter and battery efficiencies, respectively. To compensate losses associated with the above-mentioned efficiency, an optional feedback should be considered. In the following subsections, three smoothing algorithms have been described and the maximum allowable ramp value has been named r max .

Ramp Control Algorithm
The ramp control strategy analyzes the output power P g flowing to the grid. After comparing the present P g value to the one of the previous iteration increased and decreased by the maximum allowable ramp, the algorithm saturate P g between these extremes: In [5], an SoC control has been proposed in order to compensate battery losses and developed by using as feedback action in Figure 5, a proportional gain "K". Although a value of K between 2 and 8 is recommended, in our case, the system is stable with a value of K equal to 1.2. During the simulation carried out, E re f has been chosen equal to half of the battery capacity, as suggested in [5].

Moving Average Algorithm
The second analyzed algorithm smooths the power evolution by computing its mean value on a specific time window T. It could be explained by the following equation: Increasing the value of time window T increases the smoothing effect of the moving average algorithm. As in the previous strategy, in [5], a feedback node aimed to compensate the losses has been introduced. In this case, the feedback block contains the integral action reported in the following equations: In Figure 6, the block scheme of this algorithm has been shown. The feedback action allows for considering the E bat values at the beginning and at the end of the day equal to each other because it compensates losses and the mean value intrinsically satisfies this requirement. The value of T has been determined by the following equation, reported in [5]:

Step-Rate Control Algorithm
The step-rate control has been presented by Marcos et al. in [5]. This algorithm can be easily explained by the flowchart in Figure 7. As the authors mentioned in [5], this algorithm relies on the high frequency smoothing effect related to interaction near the PV power plant (e.g., 6 km). In this context, a maximum ramp of 2% in one time step dt is comparable to a maximum ramp of n · 2% in a time step equal to n · dt. Hence, the algorithm works on a specific time window that lasts n times the sampling time ∆t. It starts by evaluating if the P hcpv is growing or decreasing. In case of growing, the difference between present value and the minimum value within the time window T is compared to n times the maximum allowable ∆P. In case of decreasing, the comparison takes place considering the difference between the present value and the maximum value within T, and the same allowable ∆P. In both cases, the output needed is forced to be equal to the minimum/maximum value added/subtracted by the allowable ∆P.

Smoothing the FBM Power Output
The following section considers the application of each of the above-explained smoothing algorithm after the FBM, with the aim of running the FBM at the beginning of the day, and one of the smoothing algorithm during the day. In this case, the FBM power constraints are set equal to half of the maximum power (150 kW), while energy constraints are set equal to 30% and 70% of the SoC, instead of 10% and 90% of SoC. Smoothing algorithms have the other half of power available, while, during their application, energy constraints are again set equal to 10% and 90% of SoC. In order to better simulate the real operative conditions, the FBM works on forecasting data modified with two power drops; then, for the smoothing algorithms step, one more power drop and a 4% noise have been added to the forecast array.

Results
The presented results are related to the Ottana power plant in Sardinia. The facility consists of a Fresnel-based CSP plant (8400 m 2 of net aperture), having oil as Heat Transfer Fluid (HTF), connected to an ORC turbine (gross power of 664 kW e), a direct two-tank TES system (capacity of about 15 MWhth ), a biaxial tracking CPV power plant (429 kW e) and a Sodium-Nickel (NaNiCl) Energy Storage System (ESS) with a capacity of 430 kWh and a maximum power of 300 kW.

Forecast Results
Equations (1) and (2) for the estimation of irradiance components from the cloud cover have been verified in [12,14] on measuring stations in the US territory and using different NWP tools. The validity of the approach has also been verified for the different geographical areas and for a different weather forecast service. The method has been here applied to the cloud cover forecasts obtained for the coordinates of the Elmas (CA) airport for a reference year, with hourly resolution, and compared the estimates with the historical series of GHI radiation and DNI horizontal plane projection data from the PVGIS service (https://ec.europa.eu/jrc/en/pvgis) [17], both using the European Center for Medium-Range Weather Forecasts (ECMWF) model reanalysis database (ERA5) and the reference database derived through satellite data processing (SARAH). The optimal offset parameter value has been estimated for both equations by optimizing the Root Mean Square Error (RMSE) error. Figure 8 shows, through scatterplot, the comparison between the estimated value and the value observed for GHI and for the DNI projection for the ERA5 and SARAH datasets. Figure 9 shows the forecast errors with respect to the ERA5 database through histograms and box-plots. A good match may be noticed and also verified by the values of Mean Bias Error (MBE), Mean Absolute Error (MAE), RMSE, and correlation coefficients, reported in Table 1. The table also reports the value of the Skill Score ss of the forecast evaluated with respect to the RMSE error of a reference forecast obtained through the persistence of irradiation with a lag of 24 h, which is a good estimate for a forecast horizon of one day when no tools are available for the meteorological forecast: Only the instants for which the zenith angle is less than 90 • have been considered. The table shows the results both for the offset parameter value of the [12,14] papers and, for the optimal offset value, it has been estimated in 0.60 for GHI and 0.20 for DNI. The difference in the optimal value of the offset is mainly due to the different meteorological product used for the cloud cover estimation. The error and skill score values are in line with those found both in the works already mentioned, and, in the most recent studies on the subject, see, for example, [11,18] with variations that are related to the meteorological characteristics of the geographical areas considered.
From Figure 8, it can be noticed how a part of the points of the scatter plot appear aligned on the diagonal; these are the moments in which the weather forecast has correctly estimated a null cloud cover, so the value of the irradiation corresponds to the Clear Sky that is well approximated by the Ineichen model adopted. There is also a second marked line where the points are aligned, close to the diagonal. These are points where the estimated cloud cover is very low, but the weather model can not correctly predict Clear Sky conditions or it overestimates the Cloud Cover coefficient. The results in Figure 9 show a skewness in the distribution of prediction errors, in accordance with the results of previous works. In the same figure, the box-plots show how prediction errors clearly depend on the time of day, due to varying available Clear Sky irradiation. The month-specific box-plots highlight a more uncertain forecast in the spring months due to higher cloudiness compared, for example, to the summer months, where the forecast error is, on average, lower despite higher Clear Sky irradiation. The figure also shows how forecast errors for DNI are significantly higher compared to GHI prediction errors. An improvement in forecasting can clearly be accomplished through better upstream weather forecast output, if available, by obtaining more detailed cloud cover estimation, especially for low cloudiness. A second methodology consists in the post processing of these results, which necessarily passes from a measurement of the radiation on site, possibly carried out in real time in order to exploit the feedback information on the differences between forecast and measurement. Further improvements for very short term forecasting can be achieved through the processing of satellite or sky-camera images installed on site, using optical flow techniques. In any case, the skill score values detected already demonstrate in these conditions the benefit of weather forecasting tools for radiation prediction.

Linear Model Results
The linear regression has been carried out in the Matlab environment. The first step consists of determining the solar position by the algorithm in [16], then calculating the AM value [15] and loading the measured data. In the second step, each linear regression equation has been applied in order to determine its related coefficients. In the third step, each model's power forecast has been calculated and compared to the real power production.
The model characterized by the lower Normalized Root Mean Square Error (N-RMSE) is the Hybrid Model as shown in Table 2. Linear coefficients adopted are shown in Table 3.  Table 3. Linear coefficient of the production Hybrid model adopted in this paper.

Smoothing Algorithm Results
In the following simulations, the maximum allowable ramp value has been set equal to 2%. The three presented smoothing algorithms have been tested on two different power profiles representing a sunny day and a cloudy one. In both cases, three power drops have been added, in order to test the algorithm in worst cases, and a 4% noise in order to simulate more realistic conditions. No procedure reported in [5] that regards the battery sizing has been carried out because the storage system has been previously chosen, and it is operative at the present time in the Ottana power plant. As reported in [5], in order to improve the performance of each algorithms, the initial SoC has been set equal to 50% for RC and SR and 10% (the technical minimum value) for MA. The application of Ramp-control, Moving average control, and Step-rate control on a cloudy day has been shown on the left in Figures 10-12, respectively, while Figures 10-12 on the right show the same models in a sunny day case. The same control algorithms have been applied after the application of the two thresholds FBM. Table 4 shows the main parameters of the FIR filter chosen.
As in the previous section, the following figures highlight the effect of this procedure on a power profile related to a cloudy day and to a sunny one. Figures 13-15 on the left shows the cloudy day condition and on the right the sunny day one, for Ramp-Control, Moving Average, and Step-Rate, respectively.  Figure 13. Top: Grid power flow with RC (P g , red) and without RC (P hcpv , grey). Middle: Battery power evolution (P ess , blue) and constraints (grey). Bottom: Battery SoC evolution (blue) and constraints (grey). On a cloudy day case (a), on a sunny day (b).
(a) (b) Figure 14. Top: Grid power flow with FBM-MA (P g , red) and without FBM-MA (P hcpv , grey). Middle: Battery power evolution (P ess , blue) and constraints (grey). Bottom: Battery SoC evolution (blue) and constraints (grey). On a cloudy day case (a), on a sunny day (b).
(a) (b) Figure 15. Top: Grid power flow with FBM-SR (P g , red) and without FBM-SR (P hcpv , grey). Middle: Battery power evolution (P ess , blue) and constraints (grey). Bottom: Battery SoC evolution (blue) and constraints (grey). On a cloudy day case (a), on a sunny day (b).
With the aim of better understanding the effectiveness of the twelve analyzed cases, three indicators have been chosen. The first one is the mean value of ramps, measured both on cloudy and sunny day, associated with the power profile without any control strategy and with each one of the six studied. Referring to (14), it can be expressed as: (21) Table 5 shows the results. Moreover, another parameter chosen is the number of ramp sign variation (from positive to negative and vice-versa), as an indicator of noise in the grid power profile. It can be determined by the following condition applied to the whole P g array of each control strategy: Table 6 shows the results associated with the above-mentioned parameter. The last indicator considers the difference between the maximum and the minimum value of SoC, in order to determine how much the battery has been used during the day. It can be easily determined for each strategy by the following equation: Results related to this parameter have been shown in Table 7.

Discussion
Results shown in Section 7.2 highlight how the Hybrid Model is characterized by the lowest N-RMSE. For this reason, it could be considered more suitable for modeling the Ottana power plant production on the basis of forecast data. The difference between N-RMSE could be explained considering the AM variation range and the strong wind that characterized the power plant site. Regarding the smoothing performances (Section 7.3), in all the cases studied, the application of an FBM before the smoothing algorithm improves the battery utilization. The MA and FBM-MA cases have the lowest number of ramp variations and, in addition, are characterized by a ramp mean value lower than any other methods. In particular, FBM-MA has the lowest ramp mean value with 0.35% on the cloudy day case and 0.42% on the sunny day one. Thus, FBM-MA can be chosen as the most suitable for the Ottana power plant because it guarantees the best performance on grid-side and in terms of battery utilization. The synergy between forecasting algorithms and real-time smoothing ones will be investigated further in the future. In particular, different peak-shaving or other kinds of forecast based algorithms could be used. In the future, the possibility to integrate a super capacitor could be also investigated, in order to compensate the high frequency power fluctuation instead of the battery, preserving its life cycle. The chosen linear model and smoothing algorithm will be implemented on an NI-CompactRIO on the Ottana power plant, and they will be tested on site during a long-term campaign.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: