The advent of FPV has been driven mostly by land scarcity for projects, energy security and decarbonization targets, and a loss in PV system efficiency at high operating temperatures. FPV has demonstrated great global market potential in the recent past, with enhanced technological development in photovoltaic modules and a reduction in the levelized cost of energy (LCOE) of PV energy systems [10,11]. Appraising FPV systems and projects has been a challenge, owing to the scarcity of suitable energy simulation tools for approximating the percentage increase in yield due to the cooling effect of the water surface and the different technologies employed for floaters housing the PV modules. However, research [12] correlated different heat loss factors in W/m2K to the configuration of the floating photovoltaic structure (i.e., free-standing and small/large footprint).

This study aimed to harness the quantitative benefits of coupling FPV and onshore wind facilities with hydropower plants by relating the proximity to the existing infrastructure and grid connection, technical characteristics of the electrical network, the water-saving potential of the hydro reservoir through optimal dispatch strategies, and consequently, reducing the seasonal variations of VRES. Moreover, this study utilized the time complementarity among hydropower, floating photovoltaics and onshore wind power to alleviate the current national estimated power deficit of 810 MW. This has been attributed to reduced hydropower generation, owing to low water levels in hydro reservoirs emanating from climate change-induced droughts in the past six years [13,14]. Furthermore, this study related the integration of VRES to the hydro reservoir water-saving potential by throttling down on hydropower generation in the presence of FPV and wind power. Even though there is growing interest in floating photovoltaics, there has been no systematic appraisal of the techno-economical potential in the Zambian context and Sub-Saharan Africa (SSA). This study provides the first national-level techno-economical site assessment of onshore wind and FPV potential using a combination of validated datasets, geospatial analytical tools, site-specific wind/PV energy production models and VRES grid assessment models near existing and future hydropower plants. Furthermore, this research will help in the implementation of renewable energy technologies, such as floating photovoltaics and onshore wind power, to help increase electricity generation and supply. The study will contribute to closing the data gaps that have existed in this field of study in Zambia. To put this into perspective, the existing national grid code does not address the technical requirements (i.e., rate of frequency change, low/high voltage fault ride-through, the extent of reactive power support, etc.) of integrating VRES into the network. Therefore, this paper also addresses the nature and depth of technical studies that will have to be completed in the future to bridge this gap and thus enhance participation from independent power producers. Moreover, the research will help decision-makers to make timely and informed decisions in this area. The paper will also form a basis for further studies in the academic realm.

Therefore, the authors are highly motivated to contribute to improving the lives of all Zambians and that of neighboring citizens, by enhancing electricity access and increasing the total power generation through the adoption of renewable energy technologies such as onshore wind farms and floating photovoltaics, thus alleviating the energy poverty being faced in the region. Additionally, Zambia has the potential to enhance FOREX (foreign exchange) earnings through power exports with the interconnected SAPP countries, mitigating the chronic trade deficit with which the country has been grappling.

Against this background, the remainder of the paper is structured as follows. The subsequent section looks at the literature review; thereafter, the development of a site assessment, and the screening and ranking methodology employed in this study, are described. This ranged from site identification to the filtering and ranking of sites based on the assigned relative weight and attribute suitability scores as adopted from the literature, industry practice and stakeholder engagement. The developed methodology was then applied to a case study in Zambia. Additionally, the limitations in the site appraisal methods and tools used were highlighted. Having appraised and ranked the sites accordingly, a scoping design methodology was developed to be applied to the site with the most promising potential (i.e., highly ranked site). Furthermore, the results of the detailed case study design and formulated models were examined and discussed. Lastly, the conclusion and recommendations from the research were drawn by referencing the research outcomes, key results, study limitations and further work to be done.

#### **2. Literature Review**

#### *2.1. Overview on FPV and Onshore Wind Potential*

The evolution of FPV has, in the recent past, included the hydropower industry, owing to the opportunity for retrofitting or installing FPV panels on the abundant water

surface area of hydro dams [6,15]. To put this into perspective, hydropower represents a vital aspect of the renewable energy system and covered approximately 16.4 percent of the global electricity generation at the end of 2017, which is equivalent to 1.27 TW and 4.185 TWh of total installed capacity and generated energy, respectively, owing to the increased technological investment in the equatorial regions and China. However, the negative impact of climate change (i.e., noticeable droughts) over the past decade in some regions of the world has necessitated the rapid penetration of solar photovoltaics and wind technologies [11]. The global FPV potential for waterbodies was mapped by research with hydropower capabilities, which included electricity generation and installation capacity in terawatt-hours and gigawatts, respectively. Figure 1 below outlines the total world energy distribution in "GW" at the top and "TWh" at the bottom [16].

**Figure 1.** (**a**) Illustrating FPV capacity distribution potential in GW. (**b**) Showing electricity generation in TWh.

The percentage requirement of water-body surface area that matches the capacity of some hydropower plants in Ghana, Brazil, Malaysia, India, Turkey, Egypt, Venezuela, and Zambia is given in Table 1, which compares the various powerplants under consideration.

**Table 1.** Showing estimated reservoir size and power generation to match hydropower capacity ([1] and authors' compilation). Reproduced from [1], the name of the publisher: ESMAP 2019.

Note: \* means percentage excluding mooring (1 MW covers ~0.01 km2), including mooring (1MW covers ~0.017 km2).

To put things into perspective, and by taking Kafue Gorge Upper as an example, Table 1 shows that approximately 14% of the dam area is required to match the existing hydropower capacity of 990 MW, excluding mooring considerations, while the value increases to about 24% by including mooring.

According to research published by Rosa-Clot and Tina, Farfan and Breyer, Cazzaniga, and Nordmann et al., [11,16,17], the potential of large-scale hydro-connected photovoltaics is vastly promising, owing to photovoltaics' technological advancement, including enhanced mooring and anchoring techniques. The Longyangxia power plant in China is an example of a large-scale hydro-PV hybrid energy generation system, with a distribution of 850 MW and 1250 MW of solar PV ground-mounted and hydropower plants, respectively. This energy mix offers a time complementarity in the output by utilizing dispatchable hydropower to reduce the power variations and voltage sags of the system, due to intermittent solar power. The network energy dispatch curve is thus met by the downward or upward throttling of hydropower, depending on whether the photovoltaics output is high or low, respectively, thereby improving the reliability and enhancing the total energy generation of the system [11,18]. Scholarly analyses had brought to light the mutual benefits of FPV systems by not only reducing algae growth and evaporation but also making a reduction in the generation cost of solar PV energy, owing to the lower operating temperatures of PV panels [16,19].

With regard to wind energy reviews, the research by Local-Arantegui and Serrano-Gonzalez [20] has shown a technological evolution toward larger machines (i.e., taller towers, longer blades and high-capacity power generators). To put this into a global perspective, the size of the wind turbines in terms of hub height, rated power and rotor diameter had increased from 30 m hub height, 30 m rotor diameters and 300 kW rated power in the late 1980s, to 87.7 m hub height, 92.7 m rotor diameter and 2.1 MW rated power at the beginning of 2015. This technological evolution has been driven mainly by the process of attaining carbon neutrality, grid code integration adherence, scaling up the process to minimize reliability issues, and further cost reductions owing to the increase in capacity factor of most projects. Moreover, higher wind speeds, and consequently high energy yields, prevail at high altitudes; as such, wind turbine technology has advanced to accommodate longer heights of wind turbines (i.e., an increase in hub height and rotor diameter) [21,22]. According to the global wind energy council and Jin et al., [23], at the beginning of 2015, wind resources had become the largest and most successful renewable technology deployment, with 370 GW of global cumulative capacity. This feat was achieved in approximately 20 years. Many wind turbine configuration types have been addressed in

the literature; nonetheless, the doubly fed induction generator (DFIG) configuration stood out in terms of mainstream technological development, owing to its high energy efficiency, low power consumption and low mechanical stress [24]. The evaluation and analysis of the impact of the DFIG on system stability and reliability have become pertinent with the increase in penetration of variable renewable energy source–VRES (i.e., wind) [25–28]. Swarna et al. [29] revealed the reactive power support capability of DFIGs at the wind turbine machine terminals during instances of active power generation curtailment.

#### *2.2. Local Context Perspective*

Zambia has great solar thermal and photovoltaic application potential (i.e., 5.5 kWh/m2/day of average solar insolation, with approximately 3000 sunshine hours per annum) [30], coupled with 13 hydropower plants, accounting for 85 percent of the total installed generation capacity (2800 MW), making the nation better suited to a mix in generation sources. A recent wind resource study conducted by the World Bank revealed great wind-speed potential (i.e., from 6 to 12 m/s) in some parts of Zambia (i.e., Luangwa, Serenje, Muchinga, etc.) for utility-scale wind power generation at heights above sea level between 80 and 200 m. These heights confirm that reasonable wind speeds with the potential for wind energy occur at great altitudes. This resulted in the validation and commissioning of a wind atlas with a mesoscale resolution, based on a 2-year period of accurate wind speed measurement data taken from the 8 meteorological masts [31,32].

Moreover, Zambia aims to become a middle-income nation by the year 2030 (Vision 2030), even though the country is faced with significant challenges in the quest to achieve this feat. Some of the issues faced include a limited infrastructure for electricity evacuation, low electrification rates, and low access to clean energy technologies. With urban and rural access to electricity at 67% and 4%, respectively, translating into a national average of access to electricity of 31%, this leaves approximately 12 million people without access [33,34]. Thus, these unelectrified households depend on other fuel types for energy consumption and utilization (conventional biomass for their cooking and heating, lighting using kerosene and candles). The high dependence on biomass has resulted in huge deforestation of about 250,000–300,000 hectares per annum [35,36]. With a power consumption of approximately 706 kWh per capita, this is below expectations relative to Zambia's economic and social potential, when compared to other resource-rich countries like Namibia, Peru, South Africa and Chile, whose per capita consumption is about 2 to 3 times higher than that of Zambia. At the end of 2016, Zambia had a gross domestic product (GDP) of USD 20.5 billion for a population of about 16 million people and scooped the eighteenth rank in terms of growth economic prospects in Africa. Currently, approximately 69% of the Zambian population has a lack of access to electricity, although the 31% with access usually experience power outages, especially during the drought seasons [34]. The country's estimated 2800 MW of total installed capacity limits economic growth, mostly in drought-ridden years when the expected generation output is lower than normal. This is because eighty-five percent of the total installed capacity is hydropower, which depends on good water resource availability. The country's 3 major hydropower plants (Kafue Gorge, Kariba North Extension and Kariba North) account for 81% of electricity production. This dependency on hydropower can be ascribed to the vast water resource availability, resulting in an estimated hydropower potential of 6000 MW. However, climate change has in the recent past shrunk the dynamics of this potential by making the electrical power system susceptible to droughts. To put this into perspective, four major droughts have been experienced in Zambia in the last fifteen years, with the most recent occurring in the 2015/2016, 2016/2017 and 2019/2020 rainfall seasons. Consequently, the difficulties in load management strategies by the country's power utility companies were exacerbated in the quest to conserve water resources. This led in turn to turmoil in the national GDP, owing to reduced economic activity from the commercial, manufacturing and mining sectors [33,34,37].

This study encourages all stakeholders involved in electricity generation by promoting the use of alternative renewable energy technologies, such as onshore wind and floating photovoltaics, to enhance the capacity of electricity in the country, which is in tune with the perceived outcomes of the Zambia energy policy of 2019. Although there is a need to develop a firm and clear policy framework for effective regulation of these renewable energy technologies, this could help abate risks in project financing and enhance investor confidence. This could be key in transforming Zambia into a prosperous middle-income country by 2030, owing to the technologies' contribution in promoting sustainable and safe electricity generation for economic development and growth. The capacity and provisions to build resilient and better climate models (i.e., global circulation models) and the increase in understanding natural variability would help in enacting sound and well-informed environmental policies that tackle the existing energy challenges faced in the country and prepare for the future at the same time [34].

#### *2.3. Role of Renewable-Energy Hybrid Systems in Energy Transition (Climate Mitigation and Dispatch)*

The fight against climate change, through the attainment of carbon neutrality, has been the major motivator toward the adoption of renewable energy systems globally [18,38]. Nevertheless, concerns about system security and stability are amplified by the huge penetration of variable renewable energy sources (VRES), such as wind and solar photovoltaics, into the electrical network grids [39]. The inherent fluctuations in VRES technologies add to the uncertainty and variability in the electric power network and could negatively impact system operations if they are not addressed [40]. Li et al. [41] define uncertainty as an unanticipated change in demand and generation balance from what was forecasted, while variability is an anticipated change in the demand-generation balance. The increase in the penetration of VRES has necessitated the need to understand grid code constraints and electrical network parameters to maintain the integrity, efficiency, and reliability of the power system [42,43]. Large-scale penetration of VRES is one of the main challenges faced in modern electric power systems, owing to the complexities in the interactions between active and reactive power flows in the network, based on system design and connection characteristics, thus impacting dispatch operating costs, network losses and the voltage profile [43].

Certain scholarly analyses [18,44] found an economical operational balance between non-dispatchable (i.e., solar) and dispatchable (i.e., hydro) power sources, hence promoting the penetration of more renewable sources. Due to the benefit of increased system efficiency and enhanced energy supply balance, many countries have adopted hybrid energy systems providing a dynamic mix of two or more energy sources [45,46]. Typical hybrid energy systems include hydro–PV [47,48], hydro–wind–thermal [49], hydro–wind [50,51] and hydro–wind–PV systems [52,53]. A recent study by Maronga et al. [54] evaluated the optimal mix of PV, concentrated solar power (CSP) and storage, for a mining context in Zimbabwe. Previous research [55] found a mix of hydro and photovoltaics to be broadly used in many countries, owing to the vast spread of solar PV as a principal renewable energy source globally, and the swift regulation response of hydropower. Consequently, regions such as SSA (i.e., Zambia) that are rich in both hydropower and solar PV renewable resources are better suited in the development and deployment of hydro–PV energy systems. Regarding the dynamic and optimal mix of renewable energy sources involving solar PV, hydro and wind, research mostly focuses on resource temporal complementarity [56–58], plant operations management [18,53,58], and the optimization of system configuration [48,58–60]. A study by Beluco et al. [56] revealed a reduction in customer power outages because of the time complementarity benefits attributed to the solar PV and hydro hybrid system. Research conducted in Italy by Francois et al. [57] revealed a decrease in energy balance fluctuations, owing to the mix of solar PV and hydro (run-of-river type) power. Kougias et al. [58] were able to relate an improvement in the output of the PV–small hydro energy system by the optimization of the tilt angle and system azimuth.

Studies on hybrid energy systems involving wind, PV, and hydro aim at enhancing reliability and system flexibility by optimally dispatching the available resources. Such scholarly analyses, however, introduce errors in the modeling process by omitting to include the stochastic tendencies of solar PV and wind power [61–63]. Furthermore, by using deterministic and stochastic programming, Wei and Liu [64] tackled the uncertainties of solar PV and wind systems. The enhanced system security, coupled with limitations on the system economy and flexibility resulting from the deterministic inclusion of spinning reserve to the dispatch model, were revealed by Wei and Liu [64] and Liu et al. [53]. Dong et al. [65] and Zou [66] revealed that by adopting a structured multi-scenario perspective, a stochastic optimization problem in nature was able to be converted to a deterministic one, with the inaccuracy in the optimization output being the main trade-off.

In the recent past, the economic coordination of energy systems had employed robust and resilient optimization techniques, owing to its efficiency in excluding large-scale sampling variable datasets and probability models that have a precise distribution [67]. The random nature of VRES necessitates the adaptions between the forecasted and the actual generation of a hybrid system, so as to meet the load curve at any instant in time [52,68]. Researchers [69] had developed a method to track real-time deviations between two consecutive energy-scheduling intervals, while attaching the variability and uncertainty cost of energy. Another research study [64] incorporated the energy curtailment of solar PV and wind as a penalty cost in the scheduling.

#### **3. Methodology**

#### *3.1. Site Appraisal and Ranking Methodology*

#### 3.1.1. Overview

The decision-making process regarding the suitability and location of sites for variable renewable energy sources (solar PV and wind) utilizes geospatial parameters, mainly involving GIS models in dynamic analysis (i.e., to capture, analyze, store, manage, and manipulate spatial or geographical data) [70–73]. To aid in formulating a ranking and geospatial data interpretation methodology, such GIS modeling is usually paired with multi-criteria decision-making (MCDM) [73]. Moreover, in the late 1990s, literature in the development of VRES models started gaining traction [74,75]. Global interest in the optimal siting of solar PV and wind in the recent past has arisen, due to the quest of attaining carbon neutrality, leading to the development of generic models based on the process shown in Figure 2. Firstly, the input parameter selection, ranging from socialeconomical, technical, and environmental factors, is completed [76]. For example, ideal wind site considerations typically include the proximity to the existing electrical network (i.e., for easy grid connection), proximity to a good road network, positioning further away from protected zones (i.e., national parks or heritage land), or further away from settlements to prevent noise and flicker, good resource potential (i.e., average wind speeds and capacity factor above sea level), and being further away from the flight path to prevent interference with radar equipment near airports. Unsuitable sites are then excluded from further analysis by scoring against the model input parameters (i.e., sites with low resource potential).

**Figure 2.** Diagram showing the typical structure of the multicriteria decision method (source: [71]). Reproduced from [71], the name of the publisher: ePrints Soton 2017.

The sites that have the potential for further development and pass the filtering stage are then scored and ranking using the weighted sum method (WSM) to assess their suitability (WSM is given in the equation below) [71,73]:

$$\mathbf{A}\_{\mathbf{i}}^{\rm WSM} = \sum\_{\mathbf{j}=1}^{n} \mathbf{w}\_{\mathbf{j}} \mathbf{a}\_{\mathbf{i}\mathbf{j}} \text{ for } \mathbf{i} = 1, 2, 3, \dots \text{N} \tag{1}$$

where w is the relative parameter weighting, a is the parameter score value, and i is the attribute layer.

3.1.2. Proposed Study Methodology

The proposed study methodology for assessing site suitability was confirmed after stakeholder engagements (i.e., local experts, power utility) and extensive reviews from the literature [77–80]. The sites of interest included 5 reservoir-type, 2 pondage-type and 7 run-of-river (RoR)-type hydro plants, as shown in Table 2.

**Table 2.** Showing the identification of the hydro sites under study.

The developed methodology for the placement of wind turbines and FPV near hydropower sites is illustrated in the flowchart given in Figure 3, below.

**Figure 3.** Diagram showing a proposed methodology flowchart of the study area shown.

3.1.3. Criteria Hierarchy Structure

#### Optimal FPV Site

A two-stage approach was utilized in the selection of FPV sites, namely, screening and filtering (stage 1), and ranking and scoring (stage 2), as shown in Figure 4. The filtering stage looked at the capacity factor, distance to the grid, water surface area and distance to protected zones as the model input parameters [81–83]. The scoring and ranking stage included the relative weight (r.w.) distribution of the energy export (20% r.w.), ease of access (15% total r.w.), demand (5% total r.w.) and floating PV potential (60% total r.w.) [84–87].

#### Optimal Wind Site

The selection of onshore wind sites utilized a two-stage approach, namely, filtering and screening (stage 1) and scoring and ranking (stage 2), which is similar to the process employed under FPV (please refer to Figure A1 in Appendix A). The filtering stage looked at the distance to the grid, distance to the protected zone, wind speed, capacity factor, the security risk of installation, and noise and flicker considerations due to proximity to buildings and settlements, as the model input parameters [88–95]. The second stage (ranking and scoring) included the distribution of the energy export (20% of r.w.), ease of access (15% of total r.w.), demand (5% of total r.w.) and wind potential (60% of total r.w.).

#### 3.1.4. Site Attribute Suitability Score

Adopted from previous research, three site-attribute suitability tables were developed (shown in Appendix A); these included onshore wind, floating photovoltaics and hybrid suitability, looking at the balanced parameters of FPV and wind [72,77–80,96–105]. Since the wind potential is less pronounced than PV in Zambia, the relative weight for the wind was set lower than that of FPV in the balanced suitability ranking.

#### 3.1.5. Methodology Limitations

The weighted sum method (WSM) was applied without having insight about the assigned relative weight to the attribute layers and the layer combination procedures [73]. On the other hand, the analytic hierarchy procedure (AHP) was able to mitigate raised concerns regarding the WSM [77], although the models remain sensitive to the adopted relative weighting, as evidenced by planning permission refusal for some high-level projects within the United Kingdom. Van Rensburg et al. [106] were able to address the weighting concerns of input parameters by establishing the relationship between significant parameters influencing the quantitative assessment-based decision and the project receiving planning permission. This was then coupled with GIS modeling to assess the geospatial parameters of influence in the UK [71].

To mitigate concerns raised about the weighted sum method, the proposed study looked at a wide range of input parameters that include environmental, social, climate, economic and topographical factors to attain a more pragmatic and acceptable site appraisal (screening and ranking) process. Additionally, this was done with the help of stakeholder engagement, the solicitation of local expert opinions and an extensive literature review in the decision-making process; consequently, this contributed to the reduction in uncertainties when categorizing the attribute suitability scoring scale owing to certain assumptions that were made.

Since there is no commercial floating PV and wind project currently in Zambia, there is an element of bias in the contribution to the study from stakeholders and experts on the renewable energy generation forecast plan and agenda in line with existing policies (i.e., Vision 2030, National Energy Policy 2019, and the Seventh National Development Plan). Moreover, the authors acknowledge that the proposed appraisal method is an ongoing process, and hence is prone to some fine-tuning, as stakeholders (i.e., project developers, investors) with specific interests and viewpoints come on board.

#### *3.2. Design Scoping Methodology*

#### 3.2.1. Design Methodology Formulation

The proposed energy system at Kafue Gorge Upper will comprise hydro, onshore wind, floating photovoltaics and grid load, as given in the schematic in Figure 5. The schematic shows the existing automatic generation control (AGC), excluding VRES, and the proposed Hydro-FPV-Wind daily dispatch (HFWDD) strategy. The model assumes that all the three sources of generation under consideration are coupled to the same generation bus. Moreover, the model receives inputs from the reservoir height variation "Hr(t)", reservoir inflow "Qin(t)", hydro generation schedule "PHYg(t)", water usage/consumption "QT(t)", the hydro virtual battery from saved water "Qs(t)", grid load "PLD(t)", penstock flowrate "Qp(t)", onshore wind output "PWDg(t)" and floating photovoltaic output "PPVg(t)".

**Figure 5.** Diagram showing the schematic for the hydro-FPV-wind grid-tied system (adopted from [107]). Reproduced from [107], the name of the publisher: Elsevier 2019.

#### 3.2.2. Hydro-FPV-Wind Daily Dispatch (HFWDD) Model

The objective of the HFWDD model is to balance the seasonal load characteristic curve on the grid by optimally dispatching the three generation sources (i.e., hydro, FPV and wind). This entails developing a two-stage model that addresses the technical parameters of the electrical network for any additional generation and, thereafter, optimizes the energy system using a customized dispatch algorithm (Figure 6). Firstly, the extent of wind and FPV integration on the grid that would negatively impact the network parameters (i.e., power losses, voltage magnitude and stability) is assessed in the two-stage model. This is in line with previous research [33,42,108–114]. Secondly, seasonal hourly reservoir inflows, water consumption targets, grid demand characteristics and the total generation scenarios (wind, FPV and hydro) are incorporated into the model. Subsequently, the grid load is served by prioritizing the integration of VRES [115] as readily available, followed by a downward regulation of hydro-generation at any moment. This throttling down of hydropower is equivalent to the water-saving potential (virtual hydro battery). Nevertheless, limited reservoir capacity, coupled with a reduction in the grid demand, could present storage challenges in a wet year (which is "rarely experienced"), hence necessitating the opening of the floodgates to get rid of excess water. Similar optimization and dispatch studies of RES were conducted by [111,116–126].

**Figure 6.** Diagram showing the systematic flow of the decision level to attain optimal hydro-FPVwind daily dispatch (HFWDD).

Without putting the stochastic nature of wind and FPV power under consideration, the optimization problem is the seasonal daily dispatch on a typical day, based on the minimization of the operating cost of the existing automatic generation controller at the

hydropower plant. Additionally, owing to the perceived low operational cost of wind and FPV, the optimization problem also prioritizes the dispatch of VRES over other conventional generation sources:

$$\text{Min}\_{\mathbf{k}} \to \text{Conv}\_{\text{operator}}(\mathbf{k}) \tag{2}$$

where "k" is the dispatch scenario for the day, including hydropower plant status, and Convoperate(.) is the daily operating cost of the power plant.

Virtual storage, as indicated in Figures 5 and 6, was modeled in HomerPro in Section 4.2.4. This looks at the availability of variable renewable energy sources (FPV and wind). Based on this factor, the model calculates how much hydro would have to be ramped down. The ramping down of hydro means that less water is utilized that is then available later (i.e., at night when the sun is not shining, and also at times when the wind is calm)—the dark/calm periods will require the saved water to be utilized.

Based on the relationship between hydrogeneration and the level of the reservoir, the Qs(t) is determined. This also looks at the reservoir rule curves that must not be violated to ensure optimal operation.

However, for a very wet year, which is "rarely experienced", the storage can have limitations in terms of reservoir capacity, and thus excess water is just wasted by opening the flood gates.

#### Parameter Uncertainty of VRES

Adopted from [59,127], the wind and FPV outputs can be represented as shown below:

$$\text{For wind} \rightarrow \text{P}\_{\text{WDg},\text{t}} \in \left[ \text{P}\_{\text{WDg},\text{t}(\text{pre})} - \text{P}\_{\text{WDg},\text{t}(\text{fu})\_{\text{t}}} \text{ P}\_{\text{WDg},\text{t}(\text{pre})} + \text{P}\_{\text{WDg},\text{t}(\text{fu})} \right] \tag{3}$$

$$\text{For FPV} \rightarrow \text{P}\_{\text{PVg},t} \in \left[\text{P}\_{\text{PVg},t(\text{pre})} - \text{P}\_{\text{PVg},t(\text{flu})}, \text{P}\_{\text{PVg},t(\text{pre})} + \text{P}\_{\text{PVg},t(\text{flu})}\right] \tag{4}$$

where PWDg/PVg,t(pre) is the predicted VRES output, PWDg/PVg,t(flu) is the maximum output fluctuation, and PWDg/PVg,t is the time-dependent power output of the VRES for any given day.

#### Model Objective Function

Cost parameters are considered for the different generation stages to attain the economical and optimum dispatch scenario "k". The hydro unit's generation cost Convoperate ("Cope = in short form" = CHYg) is the first stage. The second stage (C+ ope) mostly includes the hydro unit's adjustment costs CHYgΔ, curtailment costs of FPV and wind, given as CPVg(curt) and CWDg(curt), respectively [128]. Thus, the cost minimization objective function is given as:

$$\mathbf{C\_{ope}} = \mathbf{C\_{HDg}} = \sum\_{t=1}^{T} \left( \mathbf{a} \times \mathbf{P\_{HYg,t}^2} + \mathbf{b} \times \mathbf{P\_{HYg,t}} + \mathbf{c} \right) \tag{5}$$

<sup>C</sup>+ope = CPVg(curt) + CWDg(curt) + CHYg<sup>Δ</sup> = ∑<sup>T</sup> t=1 - *y*PVg(curt) × PPVg,t − PPVg,t(inject) + *y*HYg × ΔPHYg,t + *y*WDg(curt) × PWDg,t − PWDg,t(inject) (6)

> where the hydro units' power output at time "t" is PHYg,t, PWDg,t(inject) and PPVg,t(inject) are the wind and FPV injected power into the grid at time "t", respectively, ʎHYg is the hydro units adjustment penalty price, ΔPHYg,t the power output adjustment of hydro units, ʎWDg(curt) and ʎPVg(curt) are curtailment penalty prices for wind and FPV, respectively, and "a", "b" and "c" are hydro unit cost coefficients.

#### HFWDD Model Constraints

*Hydro Constraints:*

$$\mathbf{P\_{HYg(max)}} \ge \mathbf{P\_{HYg,t}} \ge \mathbf{P\_{HYg(min)}} \tag{7}$$

$$\mathbf{Q\_{HY\_{\tilde{\rm B}}t}} = \mathbf{y^{b}}\_{\rm HY\_{\tilde{\rm B}}} + \mathbf{y^{a}}\_{\rm HY\_{\tilde{\rm B}}} \times \mathbf{P\_{HY\_{\tilde{\rm B}}t}} \tag{8}$$

$$\mathbf{V\_{flow(max)}} \ge \mathbf{Q\_{HY\_{\tilde{\mathbf{y}}}t}} \ge \mathbf{V\_{flow(min)}} \tag{9}$$

$$\mathbf{Q}\_{\rm t+1} = \mathbf{Q}\_{\rm in,t} - \mathbf{Q}\_{\rm HYg,t(curt)} - \mathbf{Q}\_{\rm HYg,t} + \mathbf{Q}\_{\rm s,t} \tag{10}$$

$$\mathbf{Q}^{\text{max}} \ge \mathbf{Q}\_{\text{s},t} \ge \mathbf{Q}^{\text{min}} \tag{11}$$

$$\mathbf{Q}\_{\mathbf{s},1} = \mathbf{Q}\_{\mathbf{s}, \text{ini}} \tag{12}$$

$$\mathbf{Q}\_{\mathbf{s},\mathbf{T}} = \mathbf{Q}\_{\mathbf{s},\mathbf{term}} \tag{13}$$

where *y*<sup>b</sup> HYg and *y*<sup>a</sup> HYg are hydro water conversion coefficients, Qs,term and Qs,ini are final and initial storage values of the reservoir, Qmax and Qmin are the upper and lower reservoir storage limit at time "t", QHYg,t is the water consumption at any time "t" of the hydro unit, PHYg,t is the power output of the hydro unit at time "t", Qin,t is the inflow of reservoir at time "t"; Qs,t is the hydro reservoir storage at time "t", QHYg,t(curt) is the curtailment of the reservoir water, and Vflow(max) and Vflow(min) are the water consumption upper and lower limits in a given period.

*Power Flow Branch Constraints:*

$$\sum\_{i=1}^{\text{Ni}} \left( \mathfrak{f}\_{\text{bi}} \times \mathcal{P}\_{\text{it}} \right) \le \mathcal{S}\_{\text{b}(\text{max})} \tag{14}$$

where Sb(max) is the branch maximum capacity, "i" is the power system node identifier, "b" is the branch identifier, Ni is the number of system network nodes in total, Pi,t is the net active power injected into the ith node. fbi is the sensitivity factor of the bth node. *Power Balance Constraints:*

$$\rm P\_{WDg,t(pre)} + P\_{HVg,t} + P\_{PVg,t(pre)} = P\_{LD,t} \tag{15}$$

where PLD,t is the grid load of the system at any given time "t". *Onshore Wind Power Constraints:*

$$P\_{\text{WDg},t} \ge P\_{\text{WDg},t(\text{inject})} \ge 0 \tag{16}$$

where PWDg,t is the variable wind generator power output at time "t". *Floating PV Power Constraints:*

$$\mathbf{P\_{PVg,t}} \ge \mathbf{P\_{PVg,t}}(\text{inject}) \ge 0 \tag{17}$$

where PPVg,t the variable FPV power output at time "t".

#### **4. Results and Discussions**

#### *4.1. Application of Appraisal and Ranking Methodology*

#### 4.1.1. Stage 1—Site Screening

The floating photovoltaics site screening process involved the definition of five criteria that include a distance to protected zones greater than or equal to 500 m, a distance to existing electrical infrastructure less than or equal to 10 km, a capacity factor (CF) greater than or equal to 14%, and a water body surface area greater than or equal to 4000 m2. Against this benchmark, Zengamina, Victoria, and Lunzua run-of-river sites were excluded on account of having a surface area <4000 m2 to accommodate a commercially and economically viable FPV project. Further, wind site filtering involved the definition of six criteria that include a distance to protected zones (i.e., national parks) greater than or equal to 500 m, the security risk (i.e., war-prone area) of installation, an average wind speed value at 150 m above ground level greater than or equal to 6 m/s, noise and flicker allowance at five times the rotor diameter (5D), a distance to electrical infrastructure less than or equal to 60 km, and a capacity factor greater than or equal to 26%. Due to the security risk zone bordering the Democratic Republic of Congo, the Zengamina wind site was excluded from the list of potential sites. This is in line with the World Bank findings on mapping security risk-prone areas for the installation of wind validation masts. Table 3 summarizes the stage 1 screening and filtering process for all the FPV and wind sites.

**Table 3.** Table showing the combined stage 1 screening outcome for both FPV and onshore wind sites.

#### 4.1.2. Stage 2—Ranking and Scoring

Three ranking and scoring tables were developed; however, only the analysis for the balanced ranking of the hybrid system is presented for simplicity. Table 4 illustrates the site scoring results for a balanced ranking using the weighted sum method (WSM). The distribution of the relative weight for the various attributes is as follows: demand at 5%, ease of access at 15%, energy export at 20%, wind potential at 25% and floating photovoltaics at 35%. Taking "FPV distance to grid" as an example under the "Energy export" attribute layer, the application of the weight sum equation 1 is presented in Figure 7. The results analysis places Kafue Gorge Upper (KGU) and Itezhi-Tezhi at second and first rank, with total attribute values of 86.9% and 90%, respectively, while the least-ranked site is Chishimba, with a total attribute combined value of 70.6%.

Figure 7 presents an example of how to apply the weighted sum method. This looks at the "Energy Export attribute layer", with a focus on the "Distance of the floating photovoltaic plant from the grid". With reference to Table 4, part 1 of Tables A1–A3 in Appendix A, there are 5 attribute layers, namely, (*i* = 1) "Wind potential", (*i* = 2) "Floating PV potential", (*i* = 3) "Energy export", (*i* = 4) "Ease of access" and (*i* = 5) "Demand". These attribute layers have the following maximum weight distribution: attribute layer (*i* = 1) → 25%, (*i* = 2) → 35%, (*i* = 3) → 20%, (*i* = 4) → 15%, (*i* = 5) → 5%. Therefore, under "select input parameters" in Figure 7, (*i* = 3) represents the energy export attribute layer, with 20% as total weight. The energy export layer is further broken down in "FPV distance to grid" given a maximum weight of 5%, "Wind distance to grid", also given 5%, and "Grid capacity availability", given 10%. Under the "weigh input parameters" FPV distance to the grid is appearing as the first layer (*j* = 1) under the energy export attribute layer. This is assigned as "w1" with reference to Equation (1). Under "score each site against parameter", if the site's FPV distance from the grid is less than or equal to 2 km, then according to Table A1 in Appendix A, the suitability score for the site will be 100%. This is assigned

as "*a*31" from Equation (1). Therefore, to get the overall site score, the product between the assignment "*w*1" under "weigh input parameters" and the assignment "*a*31" under "score each site against parameter" is calculated. Multiplying "*w*1" by "*a*31" translates into 5% × 100%, yielding a value of 5%. The process is repeated for all other attribute layers and the layers contained underneath. The total site score is the summation of the wind total, FPV total, energy export total, ease of access total and demand total.

**Table 4.** Table showing balanced scoring and ranking matrix. Analysis based on sources [70,75–78,94–103] (detailed table shown in Appendix A—Tables A2 and A3).

**Figure 7.** Diagram showing the application of the weighted sum equation under balanced scoring and ranking.

#### *4.2. Application of Design Scoping Methodology—Kafue Gorge Upper Case Study*

After appraising the potential wind and FPV sites, the stakeholder (ZESCO Ltd.) was presented with the three ranking matrices (FPV, onshore wind and balanced) of the ten potential sites to choose from. The power utility opted to adopt the balanced ranking for the hybrid energy system, with Kafue Gorge Upper (KGU) as the chosen candidate site for detailed design. Even though Itezhi-Tezhi (ITT) was ranked first over KGU, which was second under the balanced scoring, the latter was chosen over the former owing to the following factors: the presence of a data validation wind mast at KGU, the presence of debris and dead trees in the ITT reservoir, the distance to the demand center (300 km from ITT, compared to 100 km for KGU), less reliability and stability of the grid at ITT, with one 220 kV line emanating from ITT compared to three 330 kV lines from KGU to the grid.

#### 4.2.1. VRES Grid Impact Study

Using the power system analysis toolbox (PSAT), the Zambian electrical power grid was modeled at a 330 kV voltage level comprising a 27–bus system. The model for the existing network had a real and reactive power load distribution of 2383 MW and 1061.8 MWVAr, respectively. The total modeled existing generation was 2530.6 MW real power and 857.4 MVAr reactive power, comprising the following power stations: Itezhi-Tezhi via Nambala, Lunzua via Kasama, Victoria Falls via Mukuni, Maamba Collieries Limited (MCL), PV plant at Lusaka South Multi-Facility Economic Zone (LSMFEZ), Kariba North Bank and Kafue Gorge Upper. For additional generation, 200 MW of VRES was later integrated and modeled, comprising 100 MW FPV and 100 MW wind at the KGU generation bus. The actual PSAT single-line diagram model for the network is shown in Appendix B. The key summary results for the grid impact of VRES are presented below.

#### Analysis of Existing Network

Figure 8 shows the voltage violations at 11 out of the 27,330 kV buses for the existing network before the addition of compensating equipment and additional generation (VRES, in this case). According to the Zambian grid code limits, the permissible and acceptable voltage should fall in the range between 313.5 and 346.5 kV, which is a tolerance of +/−5%. The total power losses of the existing network were 1475 MW real and −204.3 MVAr reactive power, as shown in the global power summary in Figure 9.

#### Analysis after VRES Integration

The addition of network compensating equipment corrected all 11 voltage violations with the integration of VRES at the KGU bus, comprising 100 MW wind and 100 MW floating photovoltaics, improving the voltage magnitude profile even further at buses near the KGU generation bus (Figure 10). The integration of 200 MW of VRES contributed to the reduction in the voltage support requirement by the network, owing to the drop in the generated reactive power by the Kafue Gorge hydro plant from 239 to 201 MVAr. The results also showed all the line flows to be within range (less than the 700 MVA maximum line capacity). Moreover, VRES integration also reduced the network's real power losses by 5 percent (from 147 MW to 140 MW), as can be seen in Figure 11.

**Figure 8.** Figure showing the existing 330 kV bus voltage profile without VRES integration and reactive compensation at Luano, Kansanshi, Lumwana, Kitwe, Kalumbila, and Chipata West network buses.

**Figure 9.** Figure showing the global power summary for the existing network.

**Figure 10.** Showing bus voltage after adding network compensation and VRES.

**Figure 11.** Showing power summary after VRES integration.

4.2.2. Hydro Modeling Results Analysis at KGU

The Kafue Gorge Upper hydro generation was modeled in iHoga using the power plant ratings provided by the national power utility (ZESCO Ltd. located in Zambia). The iHoga model for one turbine has 4% losses in penstock, 0% daily/hourly variability and 85% total turbine efficiency. The key summary results of the model taking a typical winter (i.e., June) and summer (i.e., November) month are presented in Figure 12a–d. In June, the maximum power output of the hydropower plant was 805 MW, corresponding to a discharge rate of 227.8 m3/s and reservoir level of 974.7 m above sea level, while the minimum generation output was 697 MW, corresponding to a discharge rate and level of 197.9 m3/s and 974.9 m, respectively. In the month of November, hydro generation output ranged between 648 and 712 MW, corresponding to a discharge rate and level range of

185–201 m3/s and 974.5–974.7 m, respectively. Figure 12e shows the hydro-generation hourly time-series graphs, serving a fraction of the total national grid load for the first day of each month of January, March, June, September and November. On 1 January, KGU generated 10.69 GWh of hydro, with an evening peak of 1.99 GWh, between 5 p.m. and 9 p.m., to serve about 27% of the total national grid demand for the day. KGU generated 15.6 GWh of hydropower on 1 March, with an evening peak of 3 GWh between 5 p.m. and 9 p.m., to serve about 46% of the total national electrical grid demand for the day. On 1 June, approximately 13.8 GWh of energy was generated, with an evening peak of 2.7 GWh between 5 p.m. and 9 p.m., to serve about 34.4% of the total demand for the day. On 1 September, 14 GWh of hydro-generation was produced, with an evening peak of 2.7GWh to serve 36% of the total grid demand for the day. Further, 1 November yielded 12.6 GWh of hydropower, with an evening peak of 2.4 GWh between 5 p.m. and 9 p.m., to serve 29.6% of the grid demand for the day.

**Figure 12.** *Cont*.

**Figure 12.** (**a**) Showing 3D June hydro-generation and discharge rate vs. time. (**b**) 3D June hydro-generation and reservoir level vs. time. (**c**) 3D November hydro-generation and discharge vs. time. (**d**) Day of month time series for hydro-generation vs. time. (**e**) Daily hydro output serving a fraction of grid load.

#### 4.2.3. VRES Modeling and Results Analysis

Results analysis for the modeling and design of the 100 MWac onshore wind and 116MWdc floating photovoltaics at Kafue Gorge Upper hydropower plant is presented.

#### Floating PV

Detailed design and modeling of the FPV system were performed using the Photovoltaic System (PVSYST) software. A form factor (DC/AC ratio) of 1.16 was adopted for the project, based on industry practice for Southern Africa. The system design comprises a parallel connection of eight sub-arrays. Each sub-array comprises one hundred and twenty series-strings of seventeen PV solar modules (unit photovoltaic module rating of 285 watt-peak with 72 polycrystalline cells) connected to an inverter rated at 500 kWac, with the A/C combiner box linking twenty-five inverters in parallel for 1 sub-array. Firstly, the photovoltaic module and array characteristics were analyzed, based on the results in Figures 13 and 14. The average PV module running temperature of between 10 and 65 ◦C yielded a minimum of 60 h of operation throughout the year, with a design standard irradiation of 1 kW/m2, an operating temperature range between 10 and 70 ◦C, coupled with corresponding module efficiency between 15.8 and 11%, respectively, at a given operating temperature range. However, at all irradiation levels, an increase in efficiency was observed with a decrease in temperature. Additionally, at 1 kW/m2, the PV module power output at maximum power point (MPP) was found to be 229.1 W (20 percent decrease) and 305.3 W (7 percent increase) at the highest and lowest operating temperatures, respectively.

According to Figure 15, the annual energy yield injected into the grid from the FPV system was 214.4 GWh/year. August had the most FPV energy injection into the grid, with 21.29 GW, while January had the least, with 13.59 GWh. Moreover, August and October had the most average global horizontal irradiation of 187.6 and 193.4 kWh/m2, respectively, while January and February had the least, with 155.3 and 155.2 kWh/m2, respectively. October had the highest average ambient temperature of 25.61 ◦C, while July had the least, at 16.78 ◦C. This corresponds to a monthly average system efficiency of 11.92 and 12.53 in the hottest and coldest months, respectively.

**Figure 13.** Showing FPV Array operational temperature.

**Figure 14.** Figure showing FPV efficiency vs. irradiation curves.

**Figure 15.** Figure showing FPV energy and efficiency vs. time (months).

It must be mentioned that in the case of generation power being greater than that injected, as indicated in Equation (17), the excess can still be injected as the local and regional demands are far higher. The case study only serves a fraction of the entire national demand; therefore, any excess can still be injected for local (Zambian) consumption or to meet part of the Southern African Power Pool's demand.

System optimization of the FPV system was performed by reviewing the impact of azimuth, pitch, tilt angle and ground-cover ratio (GCR) on grid-injected energy per year. The PVSYST embedded algorithm was used in the system optimization [129]. In the PVSYST software, after the case study location coordinates are entered, the meteorological data is selected from the list of databases. Thereafter, the design and system specifications are selected. The tilt angle values from 0 to 90 degrees are selected, with a sensitivity of 1◦ intervals, at an azimuth of 0 degrees. The azimuth is then changed to 180◦ while maintaining the same tilt angle inputs [129]. The objective function is maximizing the energy injected into the grid, which gives the simulation output in GWh. The embedded algorithm carries out a parametric analysis to search for the optimal point and plots the curves accordingly for all the input steps [129]. The ground cover ratio is optimized by looking at the ratio of the active area to the ground area. In PVsyst, the "active area" is the area of one module (length x width), multiplied by the number of modules, while the "ground area" is the area occupied by the PV array. PVSYST maximizes the injected output by tracking this ratio. The closer the ratio is to unity, the lower the injected energy into the grid. Regarding the "pitch", PVSYST will maximize the energy output by increasing the pitch. However, this requires sound engineering judgment in design by factoring in the land constraints for a particular project. Figure 16a shows a 2.6% increase in grid injected energy yield (from a base value of 214.4 to 220 GWh) at a GCR of 5%, while there was a steep decrease in yield of between 80 and 100% of the GCR ratio. Figure 16b shows that the maximum annual yield, between 214 and 215 GWh, is injected into the grid for tilt angles between 10 and 20◦ for the location in question. Figure 16c reveals that more energy is injected into the grid with every step increase in pitch (i.e., a pitch of 15 m yielded more energy compared to the baseline design value of 3 m). However, for practical considerations on space constraints, the scenario regarding pitch calls for careful analysis because, for a 400% increase (from 3 m to 15 m) in pitch, only a corresponding 2.5% increase in grid injected energy was obtained. Furthermore, a negative and positive sensitivity analysis of the azimuth angle from a baseline value of 0◦ yielded a reduction in grid-injected energy. This is because the baseline value was already at the optimized azimuth angle, as illustrated in Figure 16d. Additional analysis was able to compare the energy yield and performance ratio (PR) using the PVSYST software adjustments of the albedo and heat loss factor or U–value for floating photovoltaics (0.1 albedo and U-value 31 W/m2K) and ground-mounted (0.2 and 20 W/m2K) installations, in line with other research [72]. The results show that a floating PV has better performance (PR of 83.5% and energy yield of 214.4 GWh/y) compared to a ground-mounted system (PR of 79.3% and energy yield of 204.4 GWh/y) at the same location with similar design parameters (i.e., the tilt angle, azimuth, pitch, GCR). PVSYST was also used to evaluate the economics of the floating photovoltaic system. The analysis revealed that the cost of producing 214.4 GWh/year of energy at an investment cost of GBP 0.68/Wp was GBP 0.04/kWh, excluding operation and maintenance (O&M) costs. This FPV LCOE is competitive, with a value of GBP 0.0342/kWh and GBP 0.0335/kWh, obtained by Maronga et al. [54] and RES4Africa [130], respectively, for ground-mounted PV. Homerpro gave a more conservative annual yield of 175 GWh at an LCOE of GBP 0.067/kWh (including O&M); however, this was without factoring in the water albedo and heat loss factor of the PV module's floating island. This goes to show that a reduction in annual energy yield by approximately 18.4% increases the LCOE by almost 40% for the PVSyst and Homerpro cases that are highlighted. The cost summary and energy yield distribution for FPV and wind are summarized in Table 5.

**Figure 16.** (**a**) Optimization of injected energy vs. ground cover ratio. (**b**) Energy vs. panel tilt angle. (**c**) Energy vs. pitch in meters. (**d**) Injected energy vs. azimuth angle.

**Table 5.** Table showing the cost and energy production distribution (cost breakdown sources [130,131]).

Note: The hydro model in Homerpro had a capital cost of GBP 2.8/Wp with an operation and maintenance cost of GBP 0.017/Wp.

#### Onshore Wind

Homerpro and Renewables Ninja were used in a complementary fashion in the analysis of the KGU wind-farm output. Owing to the wide coverage of its dataset, Renewables Ninja was used to simulate the output potential of each wind turbine, whose design characteristics included: 129 m hub height, 142 m rotor diameter and 4 MW power rating per turbine. Thereafter, the Renewables Ninja wind speed output was exported to Homerpro to facilitate detailed analysis, to include the practical losses imposed on a typical wind farm with 25 by 4 MW turbines (i.e., wake effects, curtailment losses etc.). With the wind farm capacity density of approximately 6.2 MW/km2, an optimistic annual energy yield of 294 GWh was registered at the wind farm excluding system losses. However, Homerpro yielded a more conservative annual energy value, with a total of 8174 h of operation. In this scenario, about 167 GWh/year of energy was produced at a competitive LCOE of about GBP 0.07/kWh, as compared to the optimistic forecasted LCOE value of GBP 0.042/kWh obtained in the recent RES4Africa study about Zambia for the 2021/2022 benchmark [130]. The higher LCOE of wind compared to FPV is due to the fact that the resource potential for solar photovoltaics is pronounced, compared to wind in the Zambian context [30–32]. As illustrated from Figure 17a–c below, the total wind energy production on 1st January was 553.87 MWh with a peak energy of 217.54 MWh between 4–8 a.m. On 1 March, the total wind energy generated was 983.18 MWh with a peak energy of 186.1MWh between 6 a.m. and 10 a.m. and 207.2 MWh between 7 p.m. and 11 p.m. The 1st of June yielded 375.17 MWh of wind energy with a peak of 132.27 MWh between 7 a.m. and 12 p.m. and 98.2 MWh between 4 p.m. and 9 p.m. On 1 September, approximately 1654.2 MWh was generated from wind with a peak of 422.2 MWh between 6 a.m. and 11 a.m. Furthermore, 1809.66 MWh of wind was generated on 1 November, with a peak of 457.88 MWh between 6 a.m. and 11 a.m. and 385.8 MWh between 2 p.m. and 7 p.m. Any excess power indicated in Equation (16) is treated in a similar manner to FPV and is injected to meet additional local or regional demand.

#### 4.2.4. Optimal Dispatch of Hydro and VRES Hybrid System Details

With a Homerpro model (Figure 18) comprising customized virtual storage, a customized hydro initially modeled in iHoga, a PVSYST-based FPV system, and onshore wind values based on Renewables Ninja wind speed data, the following analysis can be made: the Homer Matlab Dispatch was implemented using a customized dispatch algorithm and utilized 5 hydro units, which is equivalent to 700 MW, 100 MWp of FPV, 100 MWp of wind. The Matlab code was used to ascertain a customized dispatch with high VRES penetration. Homerpro calls the Matlab Dispatch at the beginning of each time step in the simulation. The Matlab Dispatch has three input variables, namely, simulation\_state, simulation\_parameters and custom\_variables. The model virtual storage is dependent on available water in the dam and the floating PV and wind potential.

**Figure 17.** (**a**) June 3D wind speed and power vs. time. (**b**) November 3D wind speed and power vs. time. (**c**) Daily seasonal wind farm output.

**Figure 18.** Graphic showing a hybrid energy system schematic in Homerpro.

The Custom Virtual Hydro Battery at Kafue Gorge has a reservoir that can store an assumed maximum capacity of 20 million cubic meters of water (0.5 m rise, assuming it is operating at a minimum elevation of 974 m above sea level), which can discharge over a 173 h (20,000,000 m3/(32 m3/s <sup>×</sup> <sup>60</sup> <sup>×</sup> 60) period at a rate of 32 m3/s. The effective head is ~382 m, and the generator efficiency is ~(85–90)%; the power and energy of the Virtual Hydro Battery system during discharging can be calculated as follows:

Discharging:

Power generated = (p) × (g) × (v) × (h) × (eff)

where (p) is the density of water with a value of 1000 kg/m3, (g) is the gravitational constant of 9.81 m/s2, (v) is the flow rate in m3/s, (h) is the head of 382 m, and (eff) is the generator efficiency value of 90%.

Power Generated = 1000 × 9.81 × 32 × 382 × 0.9 ~= 108 MW.

For a 20-million-meter cube of water at a flow rate of 32 m3/s, the water utilization duration is approximately 173 h for one turbine, based on the plant rating table shown above. However, if more turbines operate to consume the stored water, the duration would be proportional to the number of units in operation. The electrical energy generated over the 173 h is given below.

Energy generated = Power generated × hours of usage

Energy generated = 108,000 kW ×173 h = 18,684,000 kWh (~18.7 GWh) Charging:

The initial charging assumes having a wet season and thus an abundant water supply, while other charging periods of the virtual battery system involve throttling down on the hydro when there is an availability of floating photovoltaics and onshore wind. The round-trip efficiency of the virtual battery is the efficiency of the turbogenerator unit, including friction losses in the penstock (assumed to be 90% total efficiency). The maximum capacity is the maximum electrical output, divided by the nominal voltage = 18,684,000 × 1000/17,500 = ~1,067,657 amp hours, this assumes the utilization of a generation voltage of 17.5 kV for storage calculations at KGU.

#### Optimal Daily Dispatch and Reservoir Water Saving

The optimal dispatch of the hybrid energy system at KGU involved the prioritizing of FPV and wind to serve the load and the excess met by hydropower. From Figure 19a below, 1.15 GWh of VRES generation dispatch translated into a reservoir water-saving potential of 9.5% (equivalent to 1.02 GWh of generation) on 1 January. According to Figure 19b, a water-saving potential of 9.7% (equivalent to 1.52 GWh of generation) was realized with a dispatch of 1.67 GWh of VRES on 1 March. Figure 19c shows a water-saving potential of 7.2% with a dispatch of 1.14 GWh of VRES on 1 June. Both 1 September and 1 November yielded better water-saving potentials of 16.8% and 18.7% with VRES dispatch of 2.52 GWh and 2.35 GWh, respectively (Figure 19d,e). Therefore, using the customized Homer-Matlab dispatch code, 4.93 TWh of annual energy consumption was served, translating into 28 percent more demand served when compared to other default dispatch strategies embedded in Homerpro. This load was met by 166 GWh/year of wind, 175 GWh/year of floating photovoltaics, in the presence of 4.59 TWh of hydrogeneration (five out of six 140 MW hydrogenator units with a 10% reserve operating margin per unit) and at a competitive levelized cost of energy of GBP 0.055/kWh. The undispatched hydro unit presents a virtual storage potential of approximately 108 MW by a 7.4% reduction in annual hydropower generation. Moreover, the water saving potential in this study excludes the added benefit of reduced evaporation owing to the presence of retrofitted solar PV panels on the hydro reservoir.

**Figure 19.** (**a**) Optimal daily dispatch of hydro, FPV and wind on 1 January to serve a fraction of the grid demand. (**b**) Optimal daily dispatch on 1 March. (**c**) Optimal daily dispatch on 1 June. (**d**) Dispatch on 1 September. (**e**) Dispatch on 1 November.

#### **5. Conclusions**

This study presented a comprehensive assessment of integrating onshore wind and floating photovoltaics that are adjacent to future and existing hydropower sites in Zambia. All the project objectives were successfully achieved, and these included site appraisal methodology formulation to score and rank possible hydropower sites for the potential addition of onshore wind and retrofitting of floating PV, development and scoping of a case study design methodology and its application. The authors presented an application of the devised screening and ranking multicriteria-based methodology for floating PV

and onshore wind, near hydro sites. The extensive data collection for stage 1, filtered off 3 sites (Lunzua, Victoria and Zengamina), thereby presented the remaining 10 sites to the stage 2 scoring and ranking process. This ranking process was developed for three scenarios, which are the balanced hybrid, floating PV and onshore wind models. The three-level scoring and ranking procedure yielded the following results: the balanced ranking placed Itezhi-tezhi and Kafue Gorge Upper (KGU) at first- and second-rank, with total attribute values of 90% and 86.9%, respectively; FPV ranking placed Itezhi-tezhi and Kafue Gorge Upper (KGU) at first- and second-rank, with total attribute values of 95% and 92.5%, respectively; the wind ranking placed Kafue Gorge Lower (KGL) and Kafue Gorge Upper (KGU) at first- and second-rank, with total attribute values of 93.8% and 83.8%, respectively. In all the three scoring and ranking levels, the Chishimba site was ranked the least. This study presents great insight for planners and prospective investors in floating photovoltaics and onshore wind as the factors influencing the suitability of the respective sites can easily be understood.

Moreover, the authors developed a scoping design methodology to be applied at any one of the 10 potential sites. The summarized methodology for the case study application includes assessing the technical parameters of the local electrical grid for integration of variable renewable energy sources (VRES), assessing current seasonal hydro generation and grid electrical demand in a year on an hour-by-hour basis, detailed assessment and design of the VRES (floating photovoltaics and onshore wind) for the chosen case study, assessing the storage potential (implied by throttling down hydro in the presence of VRES for the reservoir type), optimizing daily energy production of the system within grid constraints and ascertaining the levelized cost of the energy of the system.

The results of the case study at Kafue Gorge Upper were promising, with VRES integration potential within grid limits of 341 GWh and 508 GWh per annum, for the conservative and optimistic case, respectively. Furthermore, it is worth noting that the floating PV is not being presented as a competitor to ground-mounted systems, but rather as a complementary technology in specific applications (i.e., retrofitting on hydro reservoirs). Along with providing such benefits as reduced evaporation and algae growth, FPV systems have lower operating temperatures and potentially reduce the costs of solar energy generation. To put this into perspective, the current study using PVSYST showed that floating photovoltaics have a better energy yield compared to a ground-mounted system, as evidenced by a 7.4%, 5.8% and 4.9% increase in energy production for the freestanding, small-footprint and large-footprint FPV configurations, respectively, at a reduced generation cost of GBP 0.04/kWh.

Therefore, floating PV and onshore wind integration could present added technoeconomic benefits by fast-tracking new capacity development with opportunities for private investments (IPPs), new opportunities for the Zambian service and manufacturing sectors, power structure decentralization, owing to the wide spread of the renewable resources in the country (i.e., solar PV and wind are more diffused) compared to localized hydropower projects (usually located near large lakes and rivers).

#### **6. Recommendations**

The following future work is recommended to add more value and traction to the project research:

is the potential to conduct detailed network analysis to cover: N-1static security assessment, network fault level analysis and protection coordination, short-term and long-term frequency response, the effects of a spinning reserve, transient and dynamic stability performance, and voltage regulation during transience. Consequently, these studies will further provide the necessary technical requirements in the national grid code for large-scale VRES integration.

**Author Contributions:** Conceptualization, K.J.N. and P.G.T.; software, K.J.N.; methodology, K.J.N. and P.G.T.; formal analysis, K.J.N.; writing—original draft preparation (based on M.Sc. dissertation), K.J.N.; writing—review and editing, A.M., P.G.T. and A.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** K.J.N. received a scholarship to study for an M.Sc. in Sustainable Engineering: Renewable Energy Systems and the Environment at the University of Strathclyde from the Commonwealth Scholarship Commission (Scholar ID: ZMCS-2019-760).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study is available on request from the corresponding authors.

**Acknowledgments:** We are grateful for the support and data obtained from industry and would like to acknowledge the contribution of the following: ZESCO Ltd. (Chief Engineer System Studies and Project Design under Distribution Development—George Muyunda and Hydrology Engineer— Kelvin Kabwe) and the Rural Electrification Authority (Chief Executive Officer—Clement Silavwe). Part of this work was conducted under the Energy Systems Research Unit at the University of Strathclyde, Glasgow.

**Conflicts of Interest:** The author declares no conflict of interest but also acknowledges a Master of Science scholarship from the Commonwealth Scholarship Commission that made it possible to conduct this study.

*Energies* **2021**, *14*, 5330

### **Appendix A**

Appraisal and Ranking Methodology.

**Table A1.** Showing the site attribute score for FPV, Wind and Hybrid.

*Energies* **2021**, *14*, 5330

**Table A2.** Showing balanced scoring and ranking matrix results.

 Showing continuation of balanced scoring and ranking matrix

results.

**Figure A1.** Graphic showing the 2-stage hierarchy structure for optimal onshore wind site selection. Abbreviations: Turbines#—number of turbines, CF—capacity factor, r.w.—relative weight, D. grid—distance to grid, G. cap.—grid capacity, L. use—land use, L. own—land ownership, D. road—distance to road, D. dem.—distance to demand center, r.w.—relative weight.

#### **Appendix B**

**Figure A2.** Layout showing the PSAT model for the existing 330 kV Zambian network.

**Figure A3.** Layout showing the PSAT FPV and wind integration model for the 330 kV Zambian network.

**Table A4.** *Cont*.

#### **References**

### *Article* **Hour-Ahead Photovoltaic Output Forecasting Using Wavelet-ANFIS**

**Chao-Rong Chen 1, Faouzi Brice Ouedraogo 2,\*, Yu-Ming Chang 1, Devita Ayu Larasati <sup>1</sup> and Shih-Wei Tan <sup>3</sup>**

**Abstract:** The operational challenge of a photovoltaic (PV) integrated system is the uncertainty (irregularity) of the future power output. The integration and correct operation can be carried out with accurate forecasting of the PV output power. A distinct artificial intelligence method was employed in the present study to forecast the PV output power and investigate the accuracy using endogenous data. Discrete wavelet transforms were used to decompose PV output power into approximate and detailed components. The decomposed PV output was fed into an adaptive neurofuzzy inference system (ANFIS) input model to forecast the short-term PV power output. Various wavelet mother functions were also investigated, including Haar, Daubechies, Coiflets, and Symlets. The proposed model performance was highly correlated to the input set and wavelet mother function. The statistical performance of the wavelet-ANFIS was found to have better efficiency compared with the ANFIS and ANN models. In addition, wavelet-ANFIS coif2 and sym4 offer the best precision among all the studied models. The result highlights that the combination of wavelet decomposition and the ANFIS model can be a helpful tool for accurate short-term PV output forecasting and yield better efficiency and performance than the conventional model.

**Keywords:** PV forecasting; ANFIS; wavelet-ANFIS; wavelet decomposition; mother wavelet function

#### **1. Introduction**

Estimating and forecasting solar power output has played an influential and critical role in integrated system management and the optimal operation of solar power in highdemand periods. Therefore, this subject is of interest both to academia and to power companies. The forecasting of upcoming events is the backbone of crisis management; as soon as this target can be accomplished, integrated system management becomes accessible [1]. Several published studies have proposed methods to forecast the future power output. Exploiting each of these proposed prediction methods generally leads to some error. Accurate prediction of PV output can provide helpful information for power management in an integrated grid [1–3].

The booming power forecasting methods proposed for solar power systems in the last decade can be divided into statistical, artificial intelligence, fuzzy inference, and hybrid methods [4]. Statistical methods, for example, auto-regressive moving average and autoregressive integrated moving average, have been used in power system prediction [5,6]. Artificial intelligence (AI) methods can be used to minimize research costs and reduce computing time as reliable alternative methods for forecasting the performance of complex systems [7]. Among others, artificial intelligence methods, multilayer perceptron neural networks [8–10], radial basis function neural networks (RBF NN) [11], physical hybrid artificial neural networks [12], recurrent neural networks [13], deep neural networks [14],

**Citation:** Chen, C.-R.; Ouedraogo, F.B.; Chang, Y.-M.; Larasati, D.A.; Tan, S.-W. Hour-Ahead Photovoltaic Output Forecasting Using Wavelet-ANFIS. *Mathematics* **2021**, *9*, 2438. https://doi.org/10.3390/ math9192438

Academic Editors: Zbigniew Leonowicz, Arsalan Najafi and Michał Jasinski

Received: 23 August 2021 Accepted: 24 September 2021 Published: 1 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

fuzzy logic (FL) [15] and ANFIS [16–18] have been developed to forecast the output power of PV sources. Artificial neural network (ANN) is among the most effective techniques to forecast solar power system output. ANN modeling of PV power forecasting has recently been studied by researchers [19,20]. Numerous studies have demonstrated the efficiency and robustness of the ANN approach for cases characterized with high nonlinearity and mapping function due to its effectiveness. ANN has a powerful learning ability even with different noise levels, the ability to simplify complex functions, wide fault tolerance, and good accomplishment for data uncertainty [19–22]. The backpropagation learning using different variants algorithms, such as the Bayesian regularization, the scaled conjugate gradient, and the Levenberg–Marquardt algorithms, with tremendous advantage as mentioned above, were used in the ANN by Basurto et al. [22].

Meanwhile, fuzzy inference systems (FIS) offer more advantages than mathematical ones. The inference operation is close to human thinking logic and can efficiently handle complex non-linear systems. [23–25]. Moreover, FIS with a backpropagation algorithm that is geared towards collecting input and output betters the effectiveness. The FIS modeling and recognition toolbox creates Takagi–Sugeno fuzzy methods from data through product space fuzzy clustering [23]. FIS has become popular in the field of engineering in the last decade, finding a considerable variety of applications—for instance, forecasting, systems pattern recognition, control theory, and power systems. Yona et al. [24] suggested a fuzzy technique to forecast hourly PV-generated power. Besides the above asset, FIS can be merged with ANN to create ANFIS [25–29]. Since Jang [25] first introduced ANFIS in the early 1990s [25], its application in power system forecasting has become rather recent. Different power system forecasting techniques have been developed to improve integrated system management [26–31]. Yaïci et al. [26] highlight that ANFIS can yield high reliability forecasting of the performance of this type of energy system. However, time series decomposition would effectively improve the ANFIS model by obtaining usable knowledge at various resolution levels to enhance the forecasting performance.

Research literature that uses wavelet and ANFIS to forecast PV-generated power, setting the model input data as PV power output, has not been found during the investigation of other proposed methods. Nevertheless, various published studies in the PV system field and other fields developed forecasting models based on the conjunction of wavelet decomposition and ANFIS. One model uses the variability reduction index, gene expression programming, wavelet transform, and ANFIS to assess the produced power for a batch of PV systems spread over one square kilometer, utilizing solar irradiance data and weather conditions [28]. Because of the irregularity and complexity of PV power, the research model [28] yielded good forecast results, but it is not easy to improve the prediction accuracy. Osorio et al. [29] proposed a short-term electricity market prices forecasting model integrating wavelet transforms and ANFIS. In addition, another research paper [30] presents a model using wavelet decomposition and ANFIS to forecast water levels. Notwithstanding, this method [30] should not include the reconstructed wavelet. Their accuracy should be calculated taking data after wavelet decomposition into consideration. Stefenon et al. [31] showed in their study that the time-series decomposition with wavelet transform improves the ANFIS robustness and accuracy and reduces training time. However, they did not propose a study of the other wavelet mother functions.

As illustrated above, in many studies, numerical weather predictions (NWP), such as temperature, relative humidity, irradiance, cloud cover, wind speed, and direction, are utilized as input to forecast the output PV power. Besides, the numerical weather data correlation with PV power is very low in the case study due to their low variation in the study area. The present study proposed endogenous data (the recorded current and past time series of the generated PV plant) as input; that is, the power generation data are directly taken, implicitly incorporating the results of NWP. Accordingly, the PV-generated power data cover the real situation on-site. In this study, two- and three-hours-previous PV output data are used to forecast solar peak time output data at different time horizons from 10 min to 60 min. The proposed idea comprises wavelet transform on the input data, training, and a testing stage using ANFIS. Compared to the existing literature, the main contributions of this paper are:

This manuscript is organized as follows: The wavelet transforms (decomposition and reconstruction) model is presented in Section 2; the overview of basic ANFIS methodology is provided in Section 3; Section 4 gives an overview of the case study. In Section 5, the forecasting results are presented with different mother wavelets, such as haar, Daubechies, Symlets, and Coiffets wavelets, on forecasting accuracy. The discussion and conclusion are presented in Sections 6 and 7, respectively.

#### **2. Wavelet Transform**

Generally, PV power output data can usually contain non-linear and dynamic features in pick and fluctuations [32], which are among the principal attributes that influence endogenous PV output power forecasting accuracy. In actual conditions, low and highfrequency signals are included in PV power data [33], which can affect the learning process of the artificial intelligence method. However, the outliers and behaviors at each frequency are more easily forecasted. Accordingly, signal decomposition techniques such as wavelet transforms can be used to forecast each frequency behavior. Consequently, the latter can improve the PV-generated forecasting accuracy. Wavelet transform is a well-documented effective "scalable" technique for time series data analysis with a higher time-frequency localization. Conventionally, discrete wavelet transforms (DWT) are generally used to increase the calculation efficiency of the prediction model. DWT of time series data can be computed as [34]:

$$DWT\_{(m,n)} = 2^{-\left(\frac{m}{2}\right)} \sum\_{t=0}^{T-1} f(t)\psi\left(\frac{t - n2^n}{2^m}\right) \tag{1}$$

where *ψ*(*t*), a continuous function, designates the wavelet mother function, *f*(*t*) is the time-series function signal, and *m* and *n* are integers used to manage the wavelet dilatation and translation, respectively. Various wavelet mother functions are available, such as Harr, Daubechies, Symlets, Coiflets, Biorthogonal, Morlet, or Mexican Hat. The most popular wavelet mother functions, including Haar, Daubechies (db3), Coiflets (coif2), and Symlets (sym4), are illustrated in Figure 1.

The Mallat technique [34] is utilized as a fast-DWT method to equilibrate the wavelength and smoothness in the current literature. This computation requires less time as it is the least complex process of calculating the DWT. The original time-series function *f*(*t*) signal can be reconstructed as the following [34]:

$$f(t) = T(t) + \sum\_{m=1}^{M} \sum\_{n=0}^{2^{M-n}-1} DWT\_{(m,n)} 2^{-\left(\frac{m}{2}\right)} \left(\frac{t - n2^n}{2^m}\right) \tag{2}$$

the integers *<sup>m</sup>* and *<sup>n</sup>* are, respectively, in the range: 1 <sup>&</sup>lt; *<sup>m</sup>* <sup>&</sup>lt; *<sup>M</sup>* and 0 <sup>&</sup>lt; *<sup>n</sup>* <sup>&</sup>lt; <sup>2</sup> *<sup>M</sup>*−*<sup>n</sup>* <sup>−</sup> 1. *T*(*t*) designates the approximation subsignal at the M level. A simple format of the signal reconstruction can be computed as [34]:

$$f(t) = \overline{T}(t) \sum\_{m=1}^{M} \mathcal{W}\_m(t) \tag{3}$$

in which *Wm*(*t*) are the details subsignals that can apprehend the interpretation of the value of small features present in the data.

**Figure 1.** Sample illustration of various wavelet mother functions: (**a**) haar, (**b**) db3, (**c**) coif2, and (**d**) sym4.

#### **3. Adaptive Neuro-Fuzzy Inference System**

ANFIS is an intelligent system that combines the strong points of the principles of fuzzy logic and neural networks into a hybrid technique. ANFIS is considered as a mere data learning algorithm that uses fuzzy logic to process data inputs into the desired output through using strongly interconnected neural network processing elements and weighted information connections to map inputs into output. This approach can simulate complex non-linear mappings, and it is suitable for an accurate short-term predictions [35–37].

As illustrated in Figure 2, the ANFIS architecture uses five layers relying on Takagi– Sugeno–Kang (TSK) rules inference system. Each layer carries out distinct functions. There are five layers: the fuzzification layer, rule reasoning layer, normalization layer, defuzzification layer, and output layer. The premise and consequent parameters are the essential parameters of ANFIS. The premise parameters {*αk*, *βk*, *γk*} are included in the designed membership functions of the fuzzification layer, which is the input layer of the ANFIS network used to generate input spaces by retrieving patterns in the input data. The consequent parameters {*ρk*, *σk*, *τk*} correspond to the parameter sets of the defuzzification layer.

The premise and the consequent parameters are optimized through training. A hybrid algorithm optimizes the parameter sets of the ANFIS forecasting system in the proposed approach.

ANFIS is explained by assuming two input variables (*z*<sup>1</sup> and *z*2) and a unique output variable (*y*ˆ). The rule illustrating the relationship associating the input, the membership function, and the output can be expressed using the if–then TSK fuzzy rules illustrated in the following conditional statements [25]:

$$\text{if } z\_1 \text{ is } D\_1 \text{ and } z\_2 \text{ is } E\_1, \text{ then } \mathcal{Y}\_1 = \rho\_1 z\_1 + \sigma\_1 z\_2 + \tau\_1 \tag{4}$$

$$\text{if } z\_1 \text{ is } D\_2 \text{ and } z\_2 \text{ is } E\_{2\prime} \text{ then } \mathcal{Y}\_2 = \rho\_2 z\_1 + \sigma\_2 z\_2 + \tau\_2 \tag{5}$$

where *D*1, *D*2, *E*<sup>1</sup> and *E*<sup>2</sup> correspond to the fuzzy sets, also called linguistic labels. The function of every layer is presented as below:

Fuzzification layer (layer 1): Each individual node is adaptive in this layer. The input variables' membership functions are mapped into fuzzy sets. The function is assigned to every node as described:

$$L\_k^1 = \mu\_{A\_k}(z\_1), \text{ for } k = 1, \ 2 \tag{6}$$

$$L\_k^1 = \mu\_{B\_{k-2}}(z\_2), \text{ for } k = 3, 4 \tag{7}$$

where *L*<sup>1</sup> *<sup>k</sup>* represents the output of the node *k*th, *μAk* and *μBk* are the membership function, there exist various membership functions. The generalized bell function can be written as:

$$\mu\_{A\_k}(z\_1) = \frac{1}{1 + \left|\frac{z\_1 - \gamma\_k}{a\_k}\right|^{2\beta\_k}} \quad 0 \le \mu\_{A\_k}(z\_1) \le 1 \tag{8}$$

where *αk*, *βk*, and *γ<sup>k</sup>* represent premise parameters used to change the shape of the membership function.

**Figure 2.** Basic ANFIS diagram.

Rule layer (layer 2): Each node in the layer is shown with a non-adaptive node, labeled as Π in Figure 2, and all the incoming signals are multiplied to compute the output:

$$L\_k^2 = w\_k = \mu\_k(z\_1).\mu\_{B\_k}(z\_2), \text{ for } k = 1, 2\tag{9}$$

where the layer 2 output is given by *L*<sup>2</sup> *<sup>k</sup>* and *wk* is weight strength for the *kth* TSK rule. Every node illustrates the weight strength of a rule.

Normalization layer (layer 3): the activity level of each rule is calculated by each node. The *k*th node is normalized as *wk*. The *kth* rule determines the weight strength, then is divided by the sum of all weight strength rules and can be computed by:

$$L\_k^3 = \overline{w}\_k = \frac{w\_k}{w\_1 + w\_2}, \text{ for } k = 1, 2 \tag{10}$$

where *L*<sup>3</sup> *<sup>k</sup>* stands for the layer 3 output and *wk* the normalized weight strength.

Defuzzification layer (layer 4): The nodes of this layer are adaptable. This layer function attributes to a single node the contribution of the *kth* rule inference to the fifth layer. *L*<sup>4</sup> *<sup>k</sup>* is computed as the defuzzification layer output using the consequent parameters and is defined as:

$$L\_k^4 = \varpi\_k \eth\_k = \varpi\_k (\rho\_k z\_1 + \sigma\_k z\_2 + \tau\_k) \tag{11}$$

Output layer (layer 5): Only one fixed node composes this layer. The whole output is computed by the sum of the arriving signals from the previous layer as:

$$\mathcal{Y} = L\_k^5 = \sum\_k \overline{w}\_k \mathcal{Y}\_k = \left(\sum\_k w\_k \mathcal{Y}\_k\right) / \left(\sum\_k w\_k\right) \tag{12}$$

#### **4. Case Study**

In power systems, for enhancing the grid performance, its stability should be guaranteed at all costs. To achieve grid stability, the appropriate scheduling of the spinning reserves and demand response is vital. With the increase in the adoption of PV as a source of energy, its intermittent nature may impact grid stability. There is a need to design forecasting models with good accuracy in forecasting the PV output power to mitigate issues with the grid stability and make PV a reliable energy source. This study used all-year solar power output data from 2017 to 2020, averaging the values for 10 min interval data. The data used for the current study were recorded at a solar power plant close to Taichung in the middle of Taiwan. The first 900 days' measurement data were employed as input training data, and the last 185 last days' data were employed to assess the trained model. The model has been designed to forecast 10 min, 30 min, and 60 min ahead each day during PV generation peak time.

The proposed forecasting algorithm consists of two parts. In the first part, the input dataset of ANFIS training data are decomposed using DWT. At that point, the function operated to decompose the training dataset is applied to the test dataset for test data decomposition. The study has also investigated the accuracy of well-known wavelet mother functions. A statistical study of different wavelet mother functions is conducted to evaluate each function's performance, including haar, db2, db3, db5, db8, coif1, coif2, coif3, coif5, sym4, sym6, and sym8.

The second part corresponds to the training and testing steps using ANFIS. Before developing the model, the optimal selection of input model numbers is essential because it can significantly reduce the computational time and cost. Two different input pattern numbers are presented in this work. The main settings of the ANFIS network include the types of input and the output membership function, the number of input and output membership functions, the number of iterations (epochs), and optimization methods such as hybrid learning.

The MATLAB toolbox is used to generate the ANFIS model for the studied data. The resulting equation of each rule is obtained by applying the linear least square assessment. Fuzzy c-means (fcm) was employed as a data clustering technique. Every data point is incumbent to a cluster and determined by membership level. The number of clusters has been set to 12, with 4 partition matrix exponents in 0.01 steps. First, the optimal dataset is carried out so that the generation of the initial FIS becomes possible. In this study, a maximum of 200 epochs has been used to achieve accurate prediction.

The training data were trained to determine the parameters of the TSK-type FIS build on the hybrid learning algorithm, which comprises the integration of the least square estimator and backpropagation gradient descent, as described by Table 1. After the training stage, the developed forecasting models were performed, and the efficiency was calculated with different evaluation criteria. The architecture of the ANFIS model developed is presented in Figure 3, with n inputs (k = 12 or k = 18), 1 output (*y*ˆ), and 12 fuzzy rules (r = 12). When k = 12, the inputs represent every 10 min PV generated data of the last 2 h; meanwhile, k = 18 employs the last 3 h of data. The 12 conditional statements are used as a rule-base. The output (*y*ˆ*i*) gives the forecasted value of PV power. Data for the same training period are used to forecast the coming 10 min, 30 min, and 60 min of each day. That is, using 2- or 3 h data to train the model (from 9 am to 11 am or 8 am to 11 am in other words t − 12, and t − 18, separately) and forecast 11:10 (t + 1), 11:30 (t + 3), and 12:00 (t+6) PV power. The training and testing flowchart are, respectively summarized in Figures 4 and 5. The asterisk (\*) in Figure 5 indicates that the test dataset wavelet transforms are acquired by applying the wavelet transform function used to decompose the training dataset.

**Table 1.** The hybrid learning methodology for ANFIS.

**Figure 3.** The proposed idea for ANFIS architecture.

**Figure 4.** Wavelet-ANFIS training flowchart.

**Figure 5.** Wavelet-ANFIS testing flowchart. The asterisk (\*) indicates that the test dataset wavelet transforms are acquired by applying the wavelet transform function used to decompose the training dataset.

The same datasets were used for an ANN model as a comparison to establish the effectiveness of the current idea. The ANN model's input pattern contained 12 or 16 inputs. The output layer was designed to forecast the power output value as described above for the ANFIS model.

#### **5. Results**

#### *5.1. Forecasting Accuracy Evaluation*

Various standard error metrics were exercised to evaluate the proposed PV output power prediction strategy. The actual and forecast sequence are represented, respectively, by . *yi* and *y*ˆ*<sup>i</sup>* with N time steps. *yi* represents the maximum recorded PV power in the test dataset. Normalize root means square error (*nRMSE* (%)), mean absolute percentage error (*MAPE* (%)), means absolute error (*MAE* (kWh)), root mean square error (*RMSE* (kWh)), and the standard deviation (*STD* (kWh) criteria are used to assess the accuracy of the different forecasting model. The metrics mentioned above are computed as below:

$$nRMSE(\%) = \sqrt{\frac{1}{N} \sum\_{i=1}^{N} \left| \frac{\dot{y}\_i - \mathcal{g}\_i}{\overline{y\_i}} \right|^2} \times 100\% \tag{13}$$

$$MAPE(\%) = \frac{1}{N} \sum\_{i=1}^{N} \left| \frac{\dot{y}\_i - \mathcal{Y}\_i}{\overline{y\_i}} \right| \times 100\% \tag{14}$$

$$MAE(\text{kWh}) = \frac{1}{N} \sum\_{i=1}^{N} |\dot{y}\_i - \mathcal{g}\_i| \tag{15}$$

$$RMSE(\text{kWh}) = \sqrt{\frac{1}{N} \sum\_{i=1}^{N} \left( \dot{y}\_i - \hat{y}\_i \right)^2} \tag{16}$$

$$STD(\text{kWh}) = \sqrt{\frac{1}{N-1} \sum\_{i=1}^{N} \left(\dot{y}\_i - \mathcal{g}\_i\right)^2} \tag{17}$$

#### *5.2. Wavelet Study Results*

The performance of the wavelet-ANFIS model for various wavelets mother functions for the input datasets is computed. Figure 6 shows the PV output of the original in blue compared with the wavelet output with different mother functions (Haar, Daubechies(db3), coiflets (coif2), and symlets (sym4)). Table 2 shows the forecasting result with different wavelet mother functions. It is discernible from Table 2 that wavelet-ANFIS models using coif2, sym4, and sym6 were more highly accurate than the other mother wavelets. The wavelet-ANFIS models with sym4 gave the highest efficiency in an overall view, indicating that wavelet decomposition utilizing the mother wavelet, sym4, can noticeably enhance ANFIS forecasting models' accuracy compared with other mother wavelets functions. It should be mentioned that the value of Table 2 is the average value computed from 30-time differently shuffled data points to have a strong comparison.

#### *5.3. Forecasting Results*

The forecasting evaluation analysis was performed by analyzing the results using the forecasting accuracy evaluation metrics mentioned in Section 5.1.

All the mother wavelets described previously have been computed, and the subsequent outcomes are presented in Table 2. It is essential to mention that after the forecasting output is obtained, the reconstitution function (2) is used to reconstruct the output. The forecasting error is computed with . *yi* , which represents the original value of the PV power output before the wavelet decomposition.

**Figure 6.** The original PV power output and different wavelet mother functions.

**Table 2.** Different wavelet-ANFIS with mother function performance.

Different forecasting methods have been completed. Firstly, wavelet-ANFIS was devised to show the effectiveness of wavelet decomposition. The decomposition significantly improves the efficiency of the forecasting method. The same forecasting cases were than conducted with ANFIS without wavelet decomposition. The statistical performance has likewise been computed with the ANN forecasting model. The ANN and ANFIS model have similar input patterns. They have inputs in the range of 12 and 18, an additional 2 hidden layers with 20 nodes at each hidden layer, and one output. The Bayesian regularization algorithm of the backpropagation method is used to update the weighting value. The maximum iteration epochs are equal to 800.

Table 3 gives the comparison of wavelet-ANFIS, ANFIS, and ANN models of the different prediction scenarios. Up to 30 scenarios are conducted for each forecasting method to randomize the initial weight. Each scenario is obtained by shuffling the available data to capture their diversity for a reliable and robust comparison. The average value of all cases RMSE is illustrated in Table 2. Table 3 summarized all evaluation criteria (nRMSE, MAPE, MAE, RMSE, and STD) of the wavelet-ANFIS, ANFIS, and ANN model.

**Table 3.** The performance of the Wavelet-ANFIS compared with ANFIS and ANN model.

#### **6. Discussion**

The model simulation results for different wavelet decomposition mother functions and two different input patterns are listed in Table 2. All the examinations are estimated by the index RMSE (kWh) defined in (13). One can observe that the mother function coif2 yielded the highest efficiency among mother wavelets in the forecasting model with 12 input patterns. Regarding the 18 input patterns, the sym4 yielded the lowest amount of forecasting errors in the overall view. Nevertheless, for the 10 min forecasting, sym6 yielded the lowest number of errors, followed by the sym4 wavelet model.

Furthermore, Table 2 also indicates that the best accuracy of 10, 30, and 60 min PV power forecasting are obtained with 18 input patterns. Meanwhile, the comparison in Table 2 demonstrates that DWT decomposition computing mother wavelets, coif2, sym4, and sym6 can considerably increase ANFIS models' effectiveness compared with the ANN one.

The comparisons between the forecasts obtained using wavelet-ANFIS, the ANFIS, and ANN are listed in Table 3. The maximums forecasting RMSE are obtained for 60 min, 2.2565 <sup>×</sup> <sup>10</sup><sup>−</sup>4, 1.0610 <sup>×</sup> <sup>10</sup><sup>−</sup>3, and 2.2924 <sup>×</sup> <sup>10</sup>−<sup>3</sup> kWh for the wavelet-ANFIS, ANFIS, ANN model, respectively. In other words, the maximum forecasting error of the ANN model was about 10.16 and 4.70 times higher than the wavelet-ANFIS and ANFIS model, separately. Such results show that the wavelet-ANFIS and ANFIS models have better performance in PV output forecasting. In addition, the MAPE, MAE, and nRMSE wavelet-ANFIS and ANFIS led to better efficiency compared with the ANN prediction model in all of the forecasting scenarios presented in this paper. The results in Table 3 are the average of 30 different shuffles, and it can be concluded that in the short-term, forecasting the PV power using the ANFIS model is more efficient than doing so using the ANN.

For more investigation, the comparison between wavelet-ANFIS and ANFIS that is presented in Table 3 indicates that the 10, 30, and 60 min for the ANFIS model with 12 input patterns forecasting RMSE are, respectively, equal to 6.3196 <sup>×</sup> <sup>10</sup><sup>−</sup>4, 8.5144 <sup>×</sup> <sup>10</sup>−4, and 1.0610 <sup>×</sup> <sup>10</sup>−<sup>3</sup> kWh, which are about 3.38, 4.04, and 4.70 times higher than the waveletmodel, separately. Meanwhile, the RMSE for the previous scenarios with 18 input patterns using the ANFIS model yielded, respectively, results about 3.90, 3.66, and 4.67 times the RMSE higher than the wavelet-ANFIS model RMSE for the exact same scenarios. From Figure 7, which illustrates the error distribution of the ANFIS and wavelet-ANFIS model, it can be observed that the 10, 30, and 60 min forecasts using wavelet-ANFIS model have more error-index computations in (18) close to zero.

$$error(\text{kWh}) = \dot{y}\_i - \mathcal{y}\_i \tag{18}$$

In summary, from Tables 2 and 3 and Figure 7, the models using wavelet components used by DWT as inputs can result in greater efficiency than the ANFIS models without wavelet transforms. It can be concluded that the forecasts using the wavelet-ANFIS models are closer to the actual PV power output values than those with the ANFIS and the ANN models. Such a result indicates that the wavelet decomposition has the capacity to enhance the effectiveness of the ANFIS forecasting PV power models.

**Figure 7.** The error distribution of wavelet-ANFIS and ANFIS with 18 input patterns: (**a**) 10 min ahead forecasting using wavelet-ANFIS model; (**b**) 10 min ahead using ANFIS model; (**c**) 30 min ahead forecasting using wavelet-ANFIS model; (**d**) 30 min ahead using ANFIS model; (**e**) 60 min ahead forecasting using wavelet-ANFIS model; (**f**) 60 min ahead using the ANFIS model.

#### **7. Conclusions**

An extensive study on applying the wavelet-ANFIS method for PV output forecasting is presented in this paper. The specific aims are to develop and evaluate a properly endogenous method for PV output forecasting in Taiwan. In this work, we compared wavelet-ANFIS, ANFIS, and ANN based on miscellaneous performance indexes, including RMSE, nRMSE, MAE, MAPE and standard deviation. The results highlight that ANFIS obtains better results than the ANN model. Furthermore, the wavelet-ANFIS model yields a better accuracy in the PV output power than the ANFIS model and the ANN model. Comparison of the wavelet-based model shows that wavelet-ANFIS-coif2 and sym4 yield better performance than all other mother functions when the last 2 and 3 h of generated PV power are used as the input of the proposed model, separately. The wavelet-ANFIS model used the last 3 h of data sym4 decomposed to yield the best overall accuracy for 10, 30, and 60 min ahead of PV power. All the same, forecasting ahead by 60 min has yielded high accuracy. The outcomes of this research demonstrate that the conjunction of DWT decomposition and the ANFIS model could meaningfully better the reliability of models used in the short-term forecasting of PV output power. The results achieved in this paper demonstrate that the connective forecast with discrete wavelet decomposition and ANFIS could be an outstanding tool for the short-term forecasting of PV output power.

**Author Contributions:** Conceptualization, C.-R.C. and F.B.O.; methodology, C.-R.C., F.B.O., Y.-M.C., D.A.L. and S.-W.T.; software, C.-R.C., F.B.O., Y.-M.C. and D.A.L.; validation, C.-R.C., F.B.O., Y.-M.C., D.A.L. and S.-W.T.; formal analysis, C.-R.C., F.B.O., Y.-M.C., D.A.L. and S.-W.T.; investigation, C.-R.C., F.B.O., Y.-M.C. and D.A.L.; resources, C.-R.C., and S.-W.T.; data curation, F.B.O. and Y.-M.C.; writing original draft preparation, C.-R.C., F.B.O. and Y.-M.C.; writing—review and editing, C.-R.C., F.B.O. and S.-W.T.; visualization C.-R.C., F.B.O., Y.-M.C. and D.A.L.; supervision, C.-R.C. and S.-W.T.; project administration, C.-R.C. and S.-W.T.; funding acquisition, C.-R.C. and S.-W.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partly supported by the University System of Taipei Joint Research Program, Project no. USTP-NTUT-NTOU-107-04.

**Institutional Review Board Statement:** The study did not involve humans or animals.

**Informed Consent Statement:** The study did not involve humans.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This study is partly supported by the University System of Taipei Joint Research Program, Project no. USTP-NTUT-NTOU-107-04.

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

#### **References**

**Chung-Geon Lee 1, La-Hoon Cho 2, Seok-Jun Kim 2, Sun-Yong Park <sup>2</sup> and Dae-Hyun Kim 3,\***

**Abstract:** The continued use of fossil fuels is contributing to severe environmental pollution and the establishment of an abnormal climate; consequently, alternative renewable energy sources are being actively investigated worldwide. Further, following global trends, numerous countermeasures aimed at improving carbon neutrality, promoting sustainable agriculture, and reducing fossil fuel dependence are being implemented in the Republic of Korea. Therefore, this study was conducted to investigate the application of renewable energies for greenhouse heating in the Republic of Korea. Three hybrid systems, numbered 1–3, were constructed using a pellet boiler, hydrothermal heat pump, and solar heat collection system, respectively. Thereafter, the heating performance, combined heat efficiency, energy consumption per temperature lifting, and energy cost per temperature lifting of the systems were compared. The combined thermal efficiency results showed no significant differences. However, in terms of energy consumption and cost, hybrid system 1 demonstrated 25.7 and 24.1% savings, respectively, compared with the other systems. Moreover, based on economic analysis via the net present value and life cycle cost analysis methods, the system reduced costs by 29.2 and 27.7%, respectively, compared with conventional fossil fuel boilers. Thus, hybrid system 1 was identified as the most economical system.

**Keywords:** hydrothermal heat pump; pellet boiler; solar heat collection; hybrid heating system; greenhouse

#### **1.  Numerous international treaties, such as the United Nations Framework Convention on Climate Change, Kyoto Protocol, and Paris Agreement, aim to reduce greenhouse gas emissions worldwide. Accordingly, most advanced countries have set targets to achieve carbon neutrality by 2050 [1–6]. The increased focus on carbon neutrality has prompted extensive research on new and renewable energy technologies. For example, to actively participate in pollution reduction and voluntarily reduce estimated carbon emissions by 37% by 2030, the Republic of Korea (ROK) announced the "2030 Greenhouse Gas Reduction Roadmap" in 2018 [7]. The roadmap included fixed targets for greenhouse gas emission reduction for each industrial sector. For the agricultural sector in particular, the target was to reduce total carbon emissions by 1%. The Energy Consumption Survey conducted in the ROK in 2017 revealed that the fossil fuel dependence of the agriculture, forestry, and fishery sectors was 97.5% [8]. Moreover, the use of fossil fuels in heating systems and in powering agricultural equipment accounted for 53.4% of all the energy consumption associated with the agricultural sector [9]. Accordingly, the

**Citation:** Lee, C.-G.; Cho, L.-H.; Kim, S.-J.; Park, S.-Y.; Kim, D.-H. Comparative Analysis of Combined Heating Systems Involving the Use of Renewable Energy for Greenhouse Heating. *Energies* **2021**, *14*, 6603. https://doi.org/10.3390/en14206603

Academic Editors: Zbigniew Leonowicz, Michał Jasinski and Arsalan Najafi

Received: 16 September 2021 Accepted: 7 October 2021 Published: 13 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

ROK is promoting the "Agricultural New and Renewable Energy Utilization Efficiency Project" to accelerate the utilization of new and renewable energy sources as well as energy-saving technologies.

In accordance with this project, heating systems using pellet boilers and geothermal heat pumps are being supplied to farms in the ROK. However, several post-installation problems are associated with these thermal energy supply systems. On the one hand, pellet boilers are difficult to maintain and have fluctuating operational costs because of variations in the cost of the wood pellets that are used as fuel [10–14]. Consequently, some farmers revert to using fossil fuel boilers. On the other hand, hydrothermal heat pumps are known to have high thermal efficiency and high heat source stability in the ROK (Republic of Korea) compared with those of air heat source pumps [15,16]. Specifically, the hydrothermal heat pumps that are supplied to farmers under the project primarily use groundwater as a heat source. However, the supply of heat, in this case, is unstable, as it depends on the amount of groundwater available [17,18]. Moreover, hydrothermal heat pumps require exorbitant initial investment costs.

Owing to these issues, farmers are hesitant to rely on alternative energy sources, and the utilization of new and renewable energy in greenhouses has remained unsuccessful.

First of all, domestic and foreign papers that applied pellet boilers to greenhouse heating include the following. Carlo Bibbiani et al. (2016) [12] judged the applicability of boilers using fossil fuels and wood fuels for greenhouse heating. In the Italian peninsula, the electric power energy load was estimated at 30 W/m<sup>2</sup> (south) to 175 W/m2 (northern), but it was estimated to be 75,362.4–1,967,796 kW/m2·yr depending on the outside temperature. The flue gas produced by the boiler contains a large amount of CO2, so recycling it has the advantage of leading to increased plant production, but biomass boilers emit more NO*x*, SO*x*, VOC, PM, and ash than fossil fuels. There are disadvantages. When using a scrubber and flue gas control device, it is possible to use a more stable wood pellet boiler, and considering this, it is highly valuable to replace with a wood biomass boiler in the case of 5–100 €/m2 depending on the greenhouse area. Kang et al. (2013) [19] designed and manufactured a wood pellet boiler to obtain basic data for practical application of a wood pellet boiler system for greenhouse heating. In order to estimate the heating efficiency according to the change in heat capacity of 75,000, 100,000, and 120,000 kcal/h, the heating efficiency test was performed by controlling the amount of air flowing into the wood pellets and burner. The thermal efficiencies of 75,000, 100,000, and 120,000 kcal/h were 80.2%, 84.2%, and 81.6%, respectively, and the highest thermal efficiency was reported at 100,000 kcal/h.

Second, there are countless papers on greenhouse heating using a heat pump, but the following studies using solar heat as an auxiliary heat source are typical for heating systems that form a heat pump and a hybrid system. Hassanien et al. (2018) [20] studied a heat pump system using a vacuum tube type solar collector as an auxiliary heat source for greenhouse heating. The internal temperature for heating was set at 14 ◦C, and it was possible to cover 62%, 40%, and 78% of the required heating load in October, March, and April, respectively. In addition, the thermal efficiency of the vacuum tube type solar collector was 0.49 and the COP of the heat pump was 4.24. However, in January, the required heating load could not be fully met, and it was judged that it was due to heat loss in the thermal storage tank. Kwon et al. (2013) [21] developed a system that improves the performance of a heat pump by selectively using solar surplus heat and external air heat in the greenhouse as heat sources, and reduces carbon dioxide fertilization costs by delaying greenhouse ventilation. Using this system, the heating performance coefficient in the internal circulation mode was about 3.35, which was improved compared to 2.46 in the night external circulation mode and 2.67 in the daytime external circulation mode. However, as the greenhouse was operated without ventilation, the light transmittance was only 62% due to excessive moisture and moisture condensation. In a horticultural environment, light transmittance is the most important factor. It was mentioned that the part showing the effect of decreasing light transmittance is considered to have room for

improvement. As such, research on a heat pump system using a pellet boiler or solar heat as an auxiliary heat source has been conducted in Korea and around the world. However, it was limited to judging the suitability or application of greenhouse heating using renewable energy. Therefore, this study was intended to study the continuous use of greenhouse heating using renewable energy or a new hybrid system.

The aforementioned solar thermal system acts as an excellent heat source under moderate climatic conditions [22]; however, its application in agriculture is challenging because of the limited land area that can be used for agriculture in the ROK. In addition, the economic feasibility of applying the solar heat collection system to agriculture without government subsidies in Korea is limited [23]. Additionally, studies on the application of solar heat in agriculture in the ROK are scarce; information regarding its performance and economic properties as a renewable energy source in agriculture is lacking.

Therefore, to ensure the use of renewable energy for heating greenhouses in ROK, this study was conducted to construct hybrid systems using new and renewable energies and to test the applicability of these systems in agriculture in the ROK. Further, based on multilateral comparative analyses and economic analyses of greenhouse heating, the most suitable system was identified.

#### **2. Materials and Methods**

#### *2.1. Greenhouse Design*

A glass greenhouse with an area of approximately 90 m2 at Kangwon University in Chuncheon-si, Gangwon-do, Korea was selected in this study. This greenhouse, which is equipped with an insulating curtain as well as ventilation facilities, features a double structure with an additional internal greenhouse. The experiments were conducted with the internal greenhouse closed. The floor area of the internal greenhouse and the covering area was 68.37 and 121.44 m2, respectively. The heating load of the experimental system was calculated based on these areas, and the capacity of each piece of equipment used in the heating system was also selected in consideration of these areas (Figure 1).

(**a**) **Figure 1.** *Cont*.

**Figure 1.** Greenhouse overview. (**a**) Photo of experimental greenhouse, and (**b**) Heating device installation overview.

#### *2.2. Thermal Energy Supply Systems*

The thermal energy supply systems tested in the greenhouse included a 20,000 kcal/hclass pellet boiler (KN-23D, Kyuwon Tech, Gyeongsan, Korea), a 3RT-class hydrothermal heat pump (3RT, Inergy Technologies, Gwangsan-gu, Korea), and a 4.04 m<sup>2</sup> solar heat collecting plate with a heat collection capacity of 2230 kcal/m2·day. Three hybrid systems were designed using these systems. The thermal energy supply device comprised a fan coil unit and a tube rail. Hybrid system 1 consisted of a hydrothermal heat pump system with a pellet boiler as the heat source. In hybrid system 2, a pellet boiler was used as the main heat source, and a solar collector was used as an auxiliary heat source for the heat storage tank. Hybrid system 3 consisted of a heat pump system that included all the heating systems, with a pellet boiler and a solar panel as the main heat source and the auxiliary heat source, respectively. The combined heating system is as shown in Figure 2.

**Figure 2.** Schematic diagram of the combined heating system.

#### *2.3. Data Measurement*

In this study, a thermocouple (GTPK-02-17, GILTRON, Seoul, Korea), sensor-type flowmeter (VVX25, SIKA Dr. Siebert & Kore GmbH & Co., Kaufungen, Germany), and turbine-type flowmeter (VTH40, SIKA Dr. Siebert & KGmb GmbH & Co., Baden-Wüttemberg, Germany) were used for data collection. The data were recorded using a data logger (GL840, GRAPHTEC Co., Tokyo, Japan), and the heat energy transfer amount was calculated by measuring the temperature and flow rate in each closed-loop system. A loadcell scale (HPS-300A, CAS Co., Seoul, Korea), an integrated watt-hour meter (LD3410DRM-080, LS ELECTRIC Co., Seoul, Korea), and a solar radiation meter (Li-200R LI-COR, Inc., Lincoln, NE, USA) were used to determine the amount of input energy. For the pellet boiler, hydrothermal heat pump, and solar heating system, the input energy was the consumption of the pellets as fuel, power consumption by the compressor, and the amount of collected solar heat, respectively. The amount of input energy was measured before heating started at 17:00, which was the standard time when the experiments commenced each day. In the experimental groups involving solar heat, the input energy was measured every day from 15:00 to noon the next day, excluding three hours (12:00–15:00) for solar heat storage. To calculate the coefficient of performance (COP) of the heat pump, the amount of instantaneous power was measured using a power analyzer (DW-6092, LUTRON ELECTRONIC ENTERPRISE C, Taipei, Taiwan).

#### *2.4. Auto Control System*

The automatic control system used in this study was the Farmos program (JINONG Co., Ltd., An-yang, Korea). This automatic control system allows communication between mobile devices and PCs, making it possible to check the operation status of various actuators, including the main heating pump, and to set manual and automatic operations. The control logic was configured such that the main heating pump and heat source supply pump could detect the temperature of the heat storage tank and the heat source tank, and accordingly set the lower and upper temperature ranges required to attain the required temperature range.

#### *2.5. Overview of Experiments*

Experiments using hybrid system 1 were conducted from 2 to 4 March 2020, and experiments using hybrid systems 2 and 3 were conducted from 12 to 17 November 2020. The heating water temperature was 55 ◦C, which is the maximum discharge temperature of

the heat pump that was used in the experiment. In hybrid systems 1 and 3, the heat source temperature was set to 20–25 ◦C. The reference heating temperature of the greenhouse was based on melon, which requires a nighttime growth temperature of 18–22 ◦C and is considered a high-temperature crop compared with other crops that are grown in Korea. Further, all the experiments were conducted via automatic control based on the control logic of the designed test method and were based on the on/off control of the thermal energy supply device.

#### *2.6. Experimental Methods*

#### 2.6.1. Hybrid System 1

The overview and control logic corresponding to hybrid system 1 are shown in Figures 3 and 4, respectively. Specifically, the operation of hybrid system 1, which enabled the detection of the temperature of the heat source tank, heat storage tank, and greenhouse, began after the set temperatures of the heat source tank and the heat storage tank were achieved and maintained before 17:00. The detected room temperature was used as input for the control logic, which determines the presence or absence of heating. To heat the greenhouse, energy was released from the heat storage tank; hence, its temperature decreased. The operation of the hydrothermal heat pump supplemented the heat from the heat storage tank based on the measured temperature of the heat storage tank. Once the temperature of the heat source tank reached the lower limit of the temperature range of the heat source, the pellet boiler was then utilized as a heat source and was operated by controlling it between the circulation pump and the heat source tank. Thus, the control logic was implemented during the realization of the experiment. Further, the daily experiments were concluded at the same time as that of the energy input check, which was performed at 17:00 the next day.

**Figure 3.** Schematic of hybrid system 1.

**Figure 4.** Control logic of hybrid system 1.

#### 2.6.2. Hybrid System 2

The overview and control logic of hybrid system 2 are shown in Figures 5 and 6, respectively. The experiment involving hybrid system 2 commenced at sunrise with the determination of the solar radiation. When the temperature of the circulating water inside the collector reached approximately 60–70 ◦C, the circulation pump between the heat exchanger and heat storage tank was operated to collect solar heat. This process continued until approximately 15:00. Further, the system continuously checked whether the heating temperature varied within the fixed range by measuring the temperature of the internal greenhouse. Heating was activated when the ambient temperature fell below the lower limit of the range. Continuous heating was performed up to the upper limit temperature of the set temperature range, and once it was determined that heating was not required, the system ceased heating activity. The collected solar heat served as an auxiliary heat source and supplemented the heat lost from the heat storage tank during the day. Additionally, the collected solar heat was used to increase the temperature of the heat storage tank to more than 55 ◦C, which represents the standard temperature for heating water, when the solar radiation intensity was high.

**Figure 5.** Schematic of hybrid system 2.

**Figure 6.** Control logic of hybrid system 2.

#### 2.6.3. Hybrid System 3

The overview and control logic of hybrid system 3 are shown in Figures 7 and 8, respectively. The operation of hybrid system 3 was based on the experimental method corresponding to hybrid system 1, with an additional solar thermal collection system. In this system, heat storage was initiated when more than a certain amount of solar radiation was detected after sunrise. Further, the heat storage process lasted from approximately 12:00 to 15:00, and once the temperature fell below the required value at 17:00, heating was initiated, with the stored solar heat utilized first. The subsequent methodology was similar to that corresponding to hybrid system 1.

**Figure 7.** Schematic of hybrid system 3.

**Figure 8.** Control logic of hybrid system 3.

#### *2.7. Analysis*

2.7.1. Thermal Energy Calculation

Heat transfer was calculated to obtain the thermal efficiency of the different hybrid systems. The amount of heat transferred was calculated using the temperature difference (ΔT), flow rate (*m*), and specific heat capacity of water (*Cp*) obtained from each closed-loop system. The calculation was performed according to Equation (1).

$$\mathbf{Q} = \Delta \mathbf{T} \times \boldsymbol{m} \times \mathbf{C}\_p \times \mathbf{3600},\tag{1}$$

where Q represents the total caloric energy (kcal/h), ΔT represents the temperature difference in in/out water (◦C), *m* represents mass flow (kg/s), and *Cp* represents the specific heat capacity of water (kcal/kg· ◦C).

#### 2.7.2. Determination of the Coefficient of Performance of the Heat Pump

The unit coefficient of performance (COP) was calculated by dividing the amount of heat transferred (kcal/h) from the heat pump to the heat storage tank (*Qhst*, *Q*heat storage tank) because of compressor power consumption (*Pcpc*, *P*compressor power consumption). The calculation was conducted by converting the amount of heat transferred into kilowatts-per-hour or converting the power consumption of the compressor into kilocalories-per-hour.

$$\text{COP} = \frac{Q\_{bst}}{P\_{cpc}} \tag{2}$$

#### 2.7.3. Energy Consumption and Cost Analysis

Energy consumption and energy cost are essential factors for the comparative analysis of the developed hybrid systems. Energy consumption was also used to calculate the efficiency of the combined thermal systems. The energy consumption of the individual systems was determined as pellet consumption, power consumption, and solar heat collection. Specifically, pellet consumption was calculated by multiplying the daily consumption amount (kg), which was obtained using a scale (HPS-300A, CAS Co., Ltd., Yang Ju, Korea), by the low-level heating amount (kcal/kg) of the pellet. Power consumption (kcal/h) was calculated by multiplying the amount of electricity measured using a watt-hour meter with a unit conversion factor, and solar heat collection was calculated by multiplying the solar heat collector efficiency by the solar radiation measured using an insolation meter (Equation (3)).

*Etotal* = (Pellet consumption × LHV) + (Power consumption × UCF) + (insolation × η), (3)

where *Etotal* represents total energy consumption (kcal), LHV represents the lower heating value (kcal/kg), UCF is the unit conversion factor (power to calories), and η represents solar collector efficiency.

The energy consumption cost was calculated based on pellet and power consumption. The standard price per unit energy was based on the wood pellet unit price (0.31 USD/kg) as announced by the Korea Forest Service in June 2019 and the Korea Electric Power Corporation (KEPCO) electricity bill calculation table (0.042 USD/kWh). The energy consumption cost, based on energy consumption and energy cost, was calculated according to Equation (4).

$$EP\_{total} = \text{(Pellet consumption } \times \text{PP)} \, + \, (\text{Power consumption } \times \text{EC}), \tag{4}$$

where *EPtotal* represents energy consumption cost, (USD, \$), PP represents pellet cost, (USD/kg), and EC represents electricity charge (USD/kWh).

#### 2.7.4. Combined Thermal Efficiency Analysis

The combined thermal efficiencies of the different hybrid systems were calculated to compare and identify the system with optimal thermal efficiency. The calculations were performed according to Equation (5), which considers the total energy consumption corresponding to each system (*Einput*) and the amount of heat transferred to the heat storage tank (*Eouput*).

$$
\eta\_{combined} = \frac{E\_{output}}{E\_{input}}\,\,\,\,\,\,\tag{5}
$$

where *Eoutput* represents energy transferred to the heat storage tank, and *Einput* represents the input energy (pellet consumption, power consumption, and solar heat collection).

#### *2.8. Economic Analysis*

#### 2.8.1. Heating Load Calculation

The total cost incurred by each system over 10 years can be predicted by dividing the total energy required for 10 years, which was calculated based on the heating load, by the energy cost per unit energy. Specifically, the heating load was calculated according to Equation (6), while the cover area heat flux, ventilation area heat flux, and floor area heat flux were calculated using Equations (7)–(9), respectively [24].

$$Q\mathbb{g} = A\_{\mathbb{X}} \times [qt + qv] + A\_{\mathbb{s}} \times q\mathbb{s} \times fw \tag{6}$$

$$qt = ht \times (T\_s - T\_a) \times (1 - fr) \tag{7}$$

$$qv = hv \times \ (T\_s - T\_a) \tag{8}$$

$$
\hbar q \text{s} = \hbar \text{s} \times (T\_s - T\_a) \tag{9}
$$

The different parameters in these equations are defined in Table 1.

**Table 1.** Factors required for the calculation of the heating load.

The total energy consumption could be calculated by multiplying the maximum daily heating load by the 10-year durability period of the device as shown in Equation (10); a 12 h non-sunlight period was assumed.

TEC (kcal) = DAHL × daily data excluding July and August × 12 h × 10 year, (10)

where TEC implies total energy consumption and DAHL implies daily average heating load.

2.8.2. Economic Analysis (Net Present Value)

The net present value method was used for the comparative economic analysis of the developed hybrid systems. In this regard, the cultivation area was assumed to be 1000 m2, and the factors considered in the net present value method of economic analysis included the durability life of the device, initial investment cost, interest rate, operating cost, and depreciation amount. Further, the initial investment cost was analyzed in two parts based on the total project cost (IIC) and the actual required project cost (IIC*sel f*) borne by the

farmers according to the Renewable Energy Use Efficiency Project conducted in Korea. The initial investment cost calculation method using the present value method was expressed as shown in Equation (11).

$$\text{TPW} = \text{IIC} \times \text{CRF} \times \text{DP},\tag{11}$$

where TPW implies total present worth, IIC implies initial investment costs, CRF implies capital recovery factor, and DP implies durability period.

The capital recovery factor (CRF) used in the net present value method was based on the straight-line depreciation method, and the cash flow was assumed to follow the same trend. The resulting CRF was calculated according to Equation (12), and the interest rate was calculated according to Equation (13) based on a nominal interest rate (2%).

$$\text{CRF} = \frac{i \times (1 + i)^n}{(1 + i)^n - 1} \tag{12}$$

where *i* represents the nominal interest rate and *n* represents the applicable year.

$$\dot{q} = \left[ (\mathbf{1} + \mathbf{r}) \times (\mathbf{1} + \mathbf{p}) \right] - \mathbf{1} \tag{13}$$

where *i* represents the nominal interest rate, r represents the real interest rate, and p represents the inflation rate.

Interest expenses, income tax, and annual operating expenses were calculated using the nominal interest rate, and the total expenses incurred during the durability period, including initial investment expenses, were calculated using Equation (14) [25].

$$TC\_{10yr} = 0 \\ TPW \\ + \left[ \left\{ \sum\_{n=1}^{10} DR \times (1+i)^n + \sum\_{n=1}^{10} IT \times (1+i)^n + \sum\_{n=1}^{10} AE \times (1+i)^n \right\} \times 10 \right] \tag{14}$$

where *TC*10*yr* represents total cost for 10 years, *DR* represents debt return, *IT* represents income tax, *AE* represents annual expenses, *i* represents nominal interest rate, and *n* represents the applicable year.

The cost per unit energy for each system was determined using the total cost for 10 years calculated above and the total energy required for 10 years calculated using Equation (10). Subsequently, Equation (15), which was used to calculate energy cost (EC), was derived as follows:

$$\text{EC} = \frac{T C\_{10yr}}{T \text{EC}}.\tag{15}$$

2.8.3. Economic Analysis (Life Cycle Cost)